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A general  time-dependent  quantum  mechanical  approach  to  the  interaction 
of  visible  and  UV  light  with  extended  polyatomic  systems  is  presented  in  this 
work.  The  interaction  of  light  with  a large  polyatomic  system  is  treated  as  a 
two-step  process:  a photon  absorption  excites  electronic  transitions  in  the  target 
system,  and  this  is  followed  by  a dynamical  evolution  of  the  system  on  the  excited 
potential  energy  surface.  Transition  rates  for  photon-molecule  interaction  processes 
are  described  in  terms  of  time-correlation  functions  (TCF)  of  transition  operators. 
These  are  written  as  a product  of  TCFs  of  the  target  and  of  the  light  source.  The 
time  evolutions  in  these  two  TCFs  are  carried  out  independently,  and  the  statistical 
properties  of  the  light  source  are  explicitly  incorporated  in  the  formalism.  The 
procedure  is  developed  for  resonance  absorption-emission  and  for  scattering  of 
thermal,  chaotic  and  coherent  light. 

The  time  evolution  in  a large  polyatomic  system  is  treated  within  a molecular 
TCF  approach.  For  a general  two-surface  electronic  excitation  problem,  a trans- 
formation of  these  molecular  TCFs  from  the  real  arguments  to  complex  time  is 
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described  so  that  their  computation  by  a numerical  path-integral  method  become 
feasible.  The  complex  time  dependence  in  these  TCFs  is  described  by  the  dy- 
namical evolution  of  an  amplitude  under  the  influence  of  a frequency  operator. 
A time-dependent  self-consistent  field  (TDSCF)  approximation  is  used  for  large 
polyatomic  system  to  factor  the  molecular  TCFs  into  primary  and  secondary  region 
TCFs.  The  primary  and  secondary  regions  are  modeled  by  considering  general 
primary  motion  coupled  to  many  harmonic  degrees  of  freedom  in  the  secondary 
region.  The  strong  coupling  between  the  primary  and  the  secondary  regions  has  a 
general  dependence  on  the  variables  of  the  primary  region,  whereas  it  has  linear 
and  bilinear  terms  in  the  variables  of  the  secondary  region.  The  weak  coupling 
limit  is  considered  by  dropping  the  bilinear  terms  in  the  coupling.  The  complex 
time  propagators  for  the  secondary  region  dynamics  for  both  of  these  cases  are 
constructed  analytically,  while  an  efficient  path-integral  method  is  developed  for 
the  dynamics  of  the  primary  region. 

As  applications,  the  path-integral  method  for  the  primary  region  dynamics  is 
used  for  computing  the  electronic  absorption  spectra  of  a model  two-potential 
problem,  and  the  TDSCF  approach  to  the  dynamics  of  large  systems  is  used 
for  studying  the  photodissociation  dynamics  of  CH3I  from  ground  and  excited 
vibrational  states  of  the  ground  potential  energy  surface.  The  cross-sections  for 
both  of  these  cases  are  compared  with  other  available  works  with  good  results. 
Extensions  of  the  method  to  the  other  extended  systems,  with  many  degrees 
of  freedom,  and  involving  tunneling  and  barrier  crossing  in  the  primary  region 
dynamics,  are  briefly  discussed. 
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CHAPTER  1 
INTRODUCTION 


Theories  of  light  absorption-emission,  elastic  and  Raman  scattering  are  de- 
scribed by  time-dependent  perturbation  expressions  in  terms  of  transition  dipole 
moments,  and  polarizability  tensors  [Kramers  and  Heisenberg,  25],  respectively. 
The  interaction  of  atoms  and  molecules  with  quantized  radiation  field  was  first 
introduced  by  Dirac  [Dirac,  27].  Since  then  many  developments  have  been  made 
along  the  same  lines  and  have  included  higher  order  corrections  in  these  theories  to 
increase  accuracies  [Herzberg,  50].  The  basic  approach  of  these  developments  is  to 
concentrate  on  small  atoms,  accurately  compute  all  of  their  eigenvalues  and  eigen- 
states, and  use  them  in  a summation  according  to  the  Fermi  Golden  Rule  [Messiah, 
62]  for  computing  the  total  cross-section  for  the  required  photointeraction  process. 
In  case  of  polyatomics,  using  the  Bom-Oppenheimer  approximation  [Balint-Kurti, 
75],  the  electronic  part  of  the  Schrodinger  equation  is  first  solved  for  fixed  nuclear 
positions  to  obtain  electronic  potential  energy  surfaces;  these  are  then  used  for  the 
nuclear  part  of  the  problem.  Given  electronic  potential  energy  surfaces  for  a small 
polyatomic  system,  one  can  still  obtain  all  the  low-lying  vibrational  and  rotational 
eigenstates  and  eigenvalues  to  calculate  cross-section.  However,  cross-sections  pro- 
cesses which  involve  rearrangement  in  the  molecule  due  to  the  absorption  of  light 
are  not  easily  described  by  the  method  of  summing  over  all  the  accessible  quantum 
states  of  the  system.  This  is  mainly  because  one  needs  exact  scattering  wavefunc- 
tions  of  many  channels  in  the  excited  state  of  the  system,  and  the  number  of  such 
available  channels  can  be  very  large  even  for  small  polyatomics.  Still,  the  time- 
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independent  quantum  mechanical  method,  based  on  summation  over  all  the  state  to 
state  information,  and  known  as  the  stationary  coupled  channel  (CC)  approach,  has 
been  used  successfully  to  produce  total  cross-sections  for  light-molecule  interac- 
tion processes  involving  rearrangement  in  the  target  [Balint-Kurti  and  Shapiro,  81; 
Heather  and  Light,  83].  However,  as  the  size  of  the  system  starts  to  increase,  meth- 
ods such  as  these  involve  a larger  and  larger  number  of  available  quantum  states  of 
the  system,  and  therefore  soon  become  too  cumbersome  to  be  of  any  practical  use. 

Alternatively,  efforts  have  been  made  to  change  the  approach  from  the  time- 
independent  quantum  mechanical  theories,  which  involve  summations  over  all  of 
the  stationary  states  of  the  system,  to  the  time-dependent  quantum  mechanical 
descriptions  which  are  based  on  the  dynamical  evolution  of  the  system  from  given 
initial  conditions  to  the  required  final  conditions  [Lee  and  Heller,  79;  Berens  and 
Wilson,  81;  Thirumalai  and  Berne,  83;  Mukamel  et  al.,  85;  Srivastava  and  Micha, 
87;  Coal  son,  87].  These  methods  do  not  require  summations  over  all  the  states 
of  the  system,  and  therefore  are  better  suited  for  dealing  with  large  molecular 
systems.  The  basic  feature  of  these  approaches  is  to  start  from  the  expression  for 
the  total  cross-section  in  the  energy  domain,  use  the  Fourier  transform  definition  of 
the  energy  conserving  delta  function  to  make  a transformation  to  the  time  domain, 
and  use  the  time-dependent  Schrodinger  equation  for  the  dynamical  evolution  of 
the  system  between  given  initial  and  required  final  conditions.  In  recent  years, 
many  advances  have  been  made  in  solving  the  time-dependent  Schrodinger  equation 
for  nuclear  vibrational  and  rotational  motion  with  given  initial  conditions  for  an 
arbitrary  general  potential  in  many  dimensions  [Kosloff,  88  for  a recent  review]. 
They  can  be  catagorized  into  the  following  broad  groups:  (a)  semi-classical  methods 
[Marcus,  72;  Miller,  74;  Mukamel  and  Jortner,  76;  Micha,  83;  Billing  and  Jolicard, 
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84;  Gray  and  Child,  84;  Henriksen,  85]  are  generally  based  on  running  classical 
trajectories  for  the  nuclear  motions,  whereas  the  electronic  transitions  are  treated 
quantum  mechanically;  (b)  approximate  quantum  mechanical  approaches  based 
on  wavepacket  dynamics  [Heller,  75,  81;  Agrawal  and  Raff,  82;  Mowrey  and 
Kouri  85;  Sawada  and  Menu,  86],  where  the  mean  position  and  the  spread  of  the 
wavepacket  follow  a canonical  set  of  equations  which  are  derived  from  a time- 
dependent  variational  principle  [Heller,  81];  (c)  straightforward  integration  of  the 
time-dependent  Schrodinger  equation  by  discretizing  the  phase  space  and  using  a 
Fast  Fourier  Transform  (FFT)  algorithm  [Feit,  Fleck  and  Steiger,  82;  Kosloff  and 
Kosloff,  83;  Hellsing,  Nitzan  and  Metiu,  86;  Kosloff,  88];  and  (d)  path  integral 
approaches  [Thirumalai  and  Berne,  83;  Miller,  Schwartz  and  Tromp,  83;  Behrman 
et  al„  83;  Coalson,  Freeman  and  Doll,  85]  which  provide  complex  time  propagators 
between  given  initial  and  final  positions  in  the  coordinate  space. 

The  present  developments  in  the  time-dependent  descriptions  still  have  the 
following  problems  which  must  be  addressed. 

By  restricting  the  radiation  field  part  of  the  interaction  Hamiltonian  to  sim- 
ple monochromatic  circular  wave  forms,  these  approaches  do  not  explicitly 
incorporate  the  statistical  and  dynamical  properties  of  the  light  source  in  the 
formalism. 

2.  Semiclassical  methods  for  solving  the  time-dependent  Schrodinger  equation  are 
based  on  classical  trajectories  which  do  not  sample  the  non-classical  regions 
of  the  potential  energy  surfaces. 

3.  Approximate  quantum  mechanical  techniques  based  on  wavepacket  dynamics 
require  local  quadratic  forms  for  potentials,  hence  they  exclude  extensive 
delocalization  of  wavepackets. 
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4.  The  numerical  path-integral  methods,  for  a potential  of  arbitrary  functional 
form,  are  not  useful  in  the  calculation  of  the  required  real  time  dynamics 
because  of  the  oscillatory  nature  of  the  real  time  propagator. 

5.  Straightforward  integration  of  the  time-dependent  Schrodinger  equation  by 
discretizing  the  phase  space  and  using  FFT  for  propagation  are  practical  only 
for  problems  with  few  spatial  dimensions. 

The  time-correlation  function  approach  to  thermodynamic  and  electromagnetic 
response  properties  of  extended  molecular  systems  is  well  known  in  statistical 
mechanics  [Forster,  75;  McQuarrie,  76;  Lovesey,  80].  Recently,  there  have 
also  been  efforts  to  describe  molecular  processes  occuring  in  a target  following 
its  collision  with  a projectile,  using  a collisional  time-correlation  function  (TCF) 
approach  [Micha,  77,  79,  81,  86].  Since,  as  before,  this  approach  does  not  require 
summations  over  all  the  states  of  the  target,  it  is  useful  in  dealing  with  many 
degrees  of  freedom  in  extended  molecular  targets  [Vilallonga  and  Micha,  87;  Micha 
and  Parra,  88],  In  this  work  we  present  a formalism  and  applications  of  a time- 
dependent  quantum  mechanical  approach  to  light-molecule  interactions  in  large 
systems  with  the  following  desirable  features. 

1.  We  treat  the  radiation  field  and  the  molecular  target  of  the  interaction  Hamilto- 
nian on  equal  footing  so  that,  at  any  stage  in  the  formalism,  explicit  statistical 
and  dynamical  properties  of  the  light  source  are  easily  incorporated. 

2.  Within  a mean  field  approximation,  we  separate  the  degrees  of  freedom  of  an 
extended  molecular  system  into  a primary  region  ( containing  a few  important 
degrees  of  freedom  along  which  all  significant  dynamical  events  take  place) 
and  a secondary  region  (with  all  of  the  remaining  degrees  of  freedom). 

3.  For  a general  potential  of  arbitrary  functional  form  in  the  primary  region,  we  use 
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a complex  time  path-integral  method  which  can  incorporate  important  quantum 
effects  such  as  tunneling  and  barrier  crossing  in  this  region. 

4.  The  computational  efforts  for  the  method  do  not  increase  quadratically  with  an 
increase  in  the  range  and  the  dimensionality  of  the  potential. 

In  chapter  2 we  describe  transition  rates  for  general  light-molecule  interaction 
processes  in  terms  of  the  TCF  of  the  Lippman- Schwinger  T operator.  In  the 
dipole  approximation  to  the  light  molecule  interaction,  the  T operator  is  expanded 
into  terms  containing  factors  for  the  radiation  field  and  the  molecular  target.  The 
transition  rates  for  all  first  order  processes  such  as  resonance  absorption  and 
emission  are  written  in  terms  of  the  molecular  dipole  TCF  of  the  target,  and  the 
electric  field  TCF  of  the  light  source.  For  all  second  order  processes,  such  as  multi- 
photon absorption-emission,  elastic  and  Raman  scattering,  the  transition  rates  are 
written  in  terms  of  the  TCFs  of  the  molecular  polarizability  tensors.  The  statistical 
properties  of  the  light  source  are  incorporated  in  the  formalism,  and  specific  cases 
of  thermal,  chaotic  and  coherent  light  are  discussed. 

In  chapter  3 we  introduce  electronic  potential  energy  surfaces  by  using  the 
Bom-Oppenheimer  approximation  to  separate  electronic  and  nuclear  motions  in  a 
molecular  target,  and  write  a general  molecular  TCF  in  the  electronic  adiabatic 
representation.  We  describe  a transformation  from  real  to  complex  time  of  these 
TCFs,  and  discuss  some  of  their  general  properties.  The  complex  time  dependence 
in  the  TCF,  for  a general  light-molecule  interaction  process  involving  electronic 
transitions  in  the  target,  is  shown  to  be  equivalent  to  the  dynamical  evolution  of 
an  amplitude  (equal  to  the  initial  eigenstate  times  the  transition  dipole)  induced  by 
a frequency  operator  (an  operator  depends  upon  the  energy  difference  between  the 
upper  and  the  lower  potential  energy  surfaces  divided  by  h). 
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In  chapter  4 we  describe  the  time  evolution  of  the  amplitude  introduced  in  the 
previous  chapter  through  a norm  conserving  time-dependent  variational  principle 
(TDVP).  Considering  fast  and  localized  electronic  tranisitions  in  a large  molecular 
target,  we  use  a time-dependent  self-consistent  field  (TDSCF)  approximation  to 
factor  the  molecular  dipole  TCF  of  the  target  into  self-consistent  primary  and 
secondary  region  TCFs. 

In  chapter  5 we  describe  time-dependent  effective  frequency  operators  for  many 
intra-  and  intermolecular  photointeraction  processes  by  considering  an  arbitrary 
general  motion  in  the  primary  region  dynamics  coupled  linearly  and  bilinearly  to 
many  degrees  of  freedom  in  a secondary  region  dynamics  undergoing  harmonic 
motions.  We  show  a direct  relationship  between  our  TDSCF  approach  and  the 
phenomenological  Generalised  Langevin  Equation  (GLE)  [Mori,  65]  approach  to 
the  stochastic  dynamics  in  large  systems.  We  explicitly  prove  a Fluctuation- 
Dissipation  Relationship  (FDR)  between  the  fluctuation  force  and  the  memory 
kernel  in  our  TDSCF  approach  for  linearly  coupled  primary  and  secondary  region 
dynamics.  We  also  develop  the  formalism  for  many  degrees  of  freedom  in  the 
primary  region  dynamics  coupled  linearly  and  bilinearly  to  many  harmonic  degrees 
of  freedom  in  the  secondary  region  dynamics. 

In  chapter  6,  we  describe  details  of  an  efficient  path-integral  method  for  the 
primary  region  dynamics.  Using  a discretized  coordinate  space,  the  computational 
method  does  not  depend  upon  any  specific  functional  form  of  the  primary  region 
potential  and  also  avoids  problems  of  other  methods  which  scale  quadratically  with 
respect  to  the  range  and  the  dimensionality  of  the  potential. 

In  chapter  7 we  describe  applications  and  results  of  the  TCF  approach  to  the 
interaction  of  visible  and  UV  light  with  extended  molecular  systems.  We  first  de- 
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scribe  an  application  to  the  electronic  absorption  spectra  of  a model  one-dimensional 
problem  for  two-potential  energy  surfaces.  The  complex  time  computations  for  the 
TCF  in  a single  spatial  dimension,  using  the  complex  time  numerical  path-integral 
method,  are  performed  so  as  to  provide  a check  of  the  accuracy  of  the  numerical 
method  for  the  primary  region  dynamics.  We  study  the  photodissociation  dynamics 
of  a realistic  polyatomic  system  (CH3I)  in  an  electronically  excited  state  within  a 
TDSCF  approximation  by  factoring  the  TCF  for  the  process  into  self-consistent  pri- 
mary and  secondary  region  TCFs.  The  resulting  total  TCF  has  been  used  to  compute 
the  total  cross-sections  for  the  photodissociation  of  CH3I,  which  are  then  compared 
with  results  from  other  methods  with  good  results.  We  have  also  extended  the 
method  to  the  computations  of  total  cross-sections  for  many  vibrationally  excited 
states  of  CH3I  on  the  ground  potential  energy  surface. 

Finally,  in  chapter  8 we  discuss  advantages  and  distinctive  features  of  our 
approach  and  mention  possibilities  for  further  improvements  and  future  extensions 
of  our  work. 


CHAPTER  2 

COLLISIONAL  TIME-CORRELATION  FUNCTION  APPROACH 
TO  LIGHT-MOLECULE  INTERACTIONS 


We  treat  the  interaction  of  light  with  a molecular  system  as  a photon  scattering 
process  through  a collisional  TCF  approach  [Srivastava  and  Micha,  87].  Starting 
from  the  transition  rates  for  a collision  process  in  terms  of  the  Lippmann-Schwinger 
transition  operator  T [Messiah,  62]  in  section  2.1,  we  derive  the  transition  rates 
for  the  interaction  of  light  with  a molecular  system  in  terms  of  the  TCF  of  these 
T operators.  In  the  dipole  approximation  to  the  light-molecule  interaction,  the 
T operator  is  written  as  an  expansion  with  terms  containing  radiation  field  and 
molecular  system  factors.  For  first  order  processes  such  as  resonance  absorption  and 
emission,  transition  rates  are  derived  in  terms  of  molecular  dipole  TCF  of  the  target 
and  the  electric  field  dipole  TCF  of  the  light  source.  For  second  order  processes 
such  as  elastic  and  inelastic  (Raman)  scattering  due  to  polarization  fluctuations  in 
the  target,  transition  rates  are  written  in  terms  of  the  polarizability  tensor  TCFs 
of  the  molecular  system  and  the  polarizability  tensor  TCFs  of  the  light  source. 
In  section  2.2  we  treat  the  interaction  of  molecules  with  thermal  light,  which 
is  given  in  the  photon  number  representation  by  a Bose-Einstein  distribution  of 
the  initial  number  of  photons  in  a particular  mode.  Section  2.3  considers  the 
interaction  of  a molecule  with  a single  mode  laser.  The  laser  light  source  is  treated 
by  considering  the  radiation  field  to  be  in  the  coherent  state  representation.  The  first 
order  transition  rates  for  laser-molecule  interaction  processes  are  derived  when  a) 
the  laser  is  operating  below  the  threshold  (chaotic  light  with  a poisson  distribution 
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for  the  mean  number  of  photons),  and  b)  the  laser  is  operating  above  the  threshold 
(coherent  light  with  a statistical  distribution  for  the  phase  and  the  mean  number  of 
photons  given  by  the  experimental  conditions). 


2.1.  Transition  Rates  from  Time-Correlation  Functions 

We  consider  the  interaction  of  visible  and  UV  light  with  a molecular  system. 
The  treatment  is  general  so  that  the  molecular  system  can  be  as  large  as  a protein 
molecule,  a polymer  chain,  an  adsorbate  on  a surface,  or  a solid  surface.  We  denote 
with  Hm,  Hr,  and  H]  the  Hamiltonians  of  the  molecular  system,  the  radiation  field, 
and  their  interaction,  respectively.  The  total  Hamiltonian  is  expressed  as 

H = Hm  + Hr  + Hj.  (2.1) 


Eigenstates  | j)  and  | r),  corresponding  to  the  noninteracting  molecular  system 
Hamiltonian  Hm  and  the  radiation  field  Hamiltonian  Hr,  respectively,  are  given 
by 


Hu\i)  = Ej\j) 


(2.2) 


A*|r)  = E,  |r). 


We  treat  the  interaction  of  light  with  a molecular  system  as  photon  scattering, 
during  which  the  radiation  field  goes  from  an  initial  state  | r)  to  a final  state  | r/). 
Later  on  we  shall  describe  various  processes  like  absorption  or  Raman  scattering  by 
specifying  | r)  and  | r/).  Transition  operators  T+  describing  the  scattering  satisfy 
the  Lippmann-Schwinger  equation  [Messiah,  62] 


T+  = Hj  + Hj K .-T+, 

E + ie-  Hm  - Hr 


(2.3) 


where  e — > 0 is  a limit  to  be  taken  after  matrix  elements  of  these  operators  have 
been  calculated.  Rates  of  transition  R,  averaged  over  the  statistical  distribution  wj 
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of  the  initial  target  states,  wr  of  the  initial  radiation  states,  and  summed  over  final 
states  of  the  target  and  the  radiation,  | jr)  and  |r/>,  respectively,  are  expressed  as 

(2.4) 

»'  rr' 

where  E = ET  + Ej  is  the  total  energy,  and  the  Dirac  8 function  imposes  energy 
conservation.  By  using  the  Fourier  transform  of  the  8—  function 

6{E  -E')  = -^-r  e^E-E^  dt,  (2.5) 

2tt  h J_  oo 

and  the  completeness  relation 

j'r')U'r'\  = (2-6) 

)'  r> 

Eq.(2.3)  can  be  written  as  [Micha,  81] 

1 f+°° 

= dt  ((T+(t)'T+(0)))M,R  (2.7) 

where  the  time-dependence  in  the  T+  is  given  by 

f+  = exp(-^of)  f+exp(+^0f)  (2.8) 

ri  h 

and  H0  = Hm  + Hr,  is  the  total  Hamiltonian  for  the  molecular  system  and 
the  radiation  field  without  interaction.  The  double  bracket  <C  ...  ^>m,r  indicates 
the  quantal  average  over  the  initial  state  |j)  of  the  molecular  system  with  a 
statistical  distribution  wj,  and  over  the  initial  state  |r)  of  the  radiation  field  with 
a statistical  distribution  wT.  So  far  we  have  not  specified  the  interaction  or  the 
coupling  Hamiltonian  #/.  Treating  the  interaction  Hamiltonian  in  the  electric 
dipole  approximation, 


Hj  = -E.D 


(2.9) 
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where  D is  the  dipole  vector  of  the  molecular  system,  and  E is  the  electric  field 
vector  of  the  light  source.  The  transition  operator  T of  Eq.(2.2)  is  written  up  to 
the  second  order,  with  an  implicit  + sign,  as 


T = ~ XI  EriDv 


1 


E + ie  — Hm  — Hr 


-r- EcD, 


<UC 


(2.10) 


£ vC 

where  (£,  77,  Q = (x,  y,  z),  and  the  time-dependence  of  the  f operator  is  given  by 
Eq.(2.8).  Since  Hm  and  Hr  commute  with  each  other,  it  is  straightforward  to  take 
matrix  elements  of  the  first  term  in  Eq.(2.10)  and  to  separate  them  into  radiation  field 
and  molecular  system  factors.  In  the  second  term  we  use  the  following  expression 


1 1 r+o° 

— -7  = - / d\ 
a + b 7 r J_00 


a 2 + A2  b2  + A2  ’ 

valid  also  for  complex  values  of  a and  b to  write 


(2.11) 


E + ie  — Hm  — Hr 


= •»  r du 

K J-  00 


Mm 


Mr 


M2m  — fi2u2  M2r  — h2u2 


(2.12) 


where 


Mm  = EM  + ieM  - Hm 

(2.13  a) 

Mr  = ER  + ieM  - Hr 

(2.13  b) 

E = Em  + Er 

(2.13 c) 

e = eM  + eR, 

(2.13 d) 

and  we  have  replaced  A = ihu.  Using  these  expressions  in  the  second  term  of 
Eq.(2.10),  to  factor  it  into  a radiation  field  part  and  a molecular  system  part,  and 
the  result  in  Eq.(2.6),  we  write  the  transition  rates  R up  to  the  second  order  as 
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p + OO 


R = A I dt  {(D\,(t)D((0)))M  ({E\,(t)Et(0)))R 

J-0 O 


X) 

-tEZT 

f'  *?C  J~°° 

p+oo 


/+oo 
-oo 


du  I dt  ((D^,(t)x^(u,  0)))m 


((4(«)V(«.o)»j! 


• + /*+oo  r-f-oo 

^2^2  du>  I dt 

((i\,c(U',t)Et(0)))R 


V vC 

n2 '+°° 


+ ^EE/  du  f+°°  du'  f+°°  dt 
2 i y_00  y.oo 


7T 


vC  ^C'-7-00 


(2.14) 


((X'J'C'(u,^)Xi?c(u’0)))a/ 

((^'C'(u,,<)^c(u’°)))« 

where  the  polarizability  tensor  operator  for  the  molecular  system  is  defined  by 


X„c(u,0)-O„a^_"vKc 

the  polarizability  tensor  operator  for  the  radiation  field  is  defined  by 


(2.15a) 


Mr 


*,c(u,0)-£,^_fiV 


< » 


(2.15fc) 


and  the  time-dependence  in  these  operators  are  given  for  example  by 


X,c(U,()  = e-^"‘«(x.O)  e+S*"'  (2.16) 

The  expression  for  i?  has  four  terms,  each  of  them  corresponding  to  different 
photointeraction  processes.  The  first  term  describes  first  order  processes  such  as 
resonance  dipolar  absorption  and  emission  of  light.  The  second  and  the  third  terms 
are  cross-terms  the  contribution  of  which  will  depend  upon  the  statistical  properties 
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of  the  light  source  used  in  the  experiment.  For  thermal  light  in  the  photon  number 
representation,  they  become  zero  and  do  not  contribute  to  the  overall  transition 
rates.  In  the  later  sections,  we  shall,  however,  see  that  the  contribution  of  these 
terms  cannot  be  ignored  if  the  light  source  used  is  a coherent  light.  The  fourth  term 
describes  second  order  processes  such  as  elastic  and  Raman  scattering  mediated  by 
polarization  fluctuations  in  the  target. 

The  transition  rates  R expressed  in  Eq.(2.14),  as  opposed  to  the  standard 
expressions  [Dirac,  47;  Heitler,  54],  are  specially  useful  for  the  treatment  of  light- 
molecule  interactions  in  large  systems  because  they  do  not  involve  summations  over 
excited  states  of  the  molecular  system.  Furthermore,  Eq.(2.10),  which  is  used  to 
separate  in  the  transition  operator  the  radiation  field  and  molecular  system  variables, 
can  also  be  used  to  derive  a similar  separation  of  the  transition  operator  to  order 
higher  than  the  second.  Thus,  a similar  formulation  of  the  transition  rates  for  higher 
order  photointeraction  processes  can  be  systematically  derived. 

Lastly,  the  main  purpose  of  separating  the  TCF  of  the  T operator  into  the 
TCF  of  the  radiation  field  and  the  TCF  of  the  molecular  system  is  to  uncouple 
the  time  evolution  of  the  molecular  system  from  that  of  the  radiation  field.  This 
has  a twofold  advantage:  firstly,  one  can  use  time-dependent  quantum  mechanical 
techniques  to  describe  the  time  evolution  of  the  molecular  system  without  referring 
to  the  time  evolution  of  the  radiation  field,  and  secondly,  we  can  easily  incorporate 
the  statistical  properties  of  the  light  source  into  the  dynamics. 

2.2.  Radiation  Field  TCFs  For  Thermal  Light 


When  the  light  source  is  maintained  at  a temperature  T,  it  is  described  by  the 
Bose-Einstein  statistical  distribution  function  W8{u}a,  T),  of  the  number  of  photons 
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in  a mode  s with  energy  ha>s.  Then  it  is  easier  to  work  in  the  photon  number 
representation  of  the  radiation  field.  The  electric  field  operator  in  Eq.(2.14)  is 
decomposed  as 


Et(t)  = E±+)  (i)  + i£  \t)  (2.17) 

where 


and 


(2.18) 


£<-»(()=  (^+l(())*  (2.19) 

where  k is  the  photon  wavevector,  e ^ is  the  £ component  of  the  polarization  vector, 
and  a ^ the  annihilation  operator  of  a photon  of  energy  = hkc  and  polarization 
£.  In  the  photon  number  representation,  a general  radiation  field  state  with  several 
modes  s is  written  as  |{(£J6)Ar‘})  where  N3  is  the  number  of  photons  in  mode 
s;  then 


%+)  (*)  lUU,)*'})  = i E (^) ! % 

(2.20) 

Using  Eqs.(2. 17-2.20)  in  Eq.(2.14),  one  can  write  the  transition  rates  for  first  and 
second  order  processes.  In  the  photon  number  representation,  the  second  and  third 
terms  of  Eq.(2.14)  do  not  contribute,  because  they  do  not  conserve  the  number  of 
photons  in  the  initial  and  the  final  radiation  state.  Therefore  the  transition  rate  is 
written  as 
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R = Ru  + R22  (2.21) 

where  Rn  is  the  transition  rate  for  the  first  order  processes,  in  which  a single 
electric  field  operator  at  an  initial  time  t = 0,  is  correlated  to  its  adjoint  at  a later 
time  t.  Similarly,  R22  is  the  transition  rate  for  second  order  processes,  where  the 
time -dependent  polarization  tensor  operator,  <^(u,  0),  which  is  second  order  in  the 
electric  field  operators  of  Eq.(2.17),  is  correlated  to  its  adjoint  at  a later  time  t.  In 
the  same  notation  R\2  and  R2\  are  the  cross-terms,  the  contributions  of  which  are 
zero  in  the  photon  number  representation.  The  transition  rates  Rn  for  first  order 
processes  are  given  by 


where 


+oo 

i?ll  = ^ f dt  (t)  , 

tt'-OO 


(2.22) 


(2.23) 

is  the  molecular  dipole  TCF  of  the  extended  system,  and 


/{J|,(<)  = ((^r)(<)£<+)(  «)))*  + ((4°  W^’w))*  (2.24) 

is  the  electric  field  TCF  of  the  light  source.  The  first  term  in  Eq.(2.24)  corresponds 
to  resonance  absorption,  whereas  the  second  term  corresponds  to  spontaneous  and 
stimulated  emissions.  Other  cross-terms  in  Rn  do  not  contribute  because  they  do 
not  conserve  the  number  of  photons.  Using  Eqs.(2. 18-2.20)  in  the  above  equation 


we  write 
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26nV' 


+ (Na  + 1W  e-'^4] 


(2.25) 


and  the  transition  rate  R\\  is  given  by 


^ = ^E»-EE  / * w «»ko  (2-26) 

s it,  £ o 

2i?e«£>£  (()'  £>£  (0)))m  + jee{e-'"‘.‘((i>£  (<)’  i><  ((>)))«}] 

where  the  statistical  probability  distribution  function  for  a thermal  light  Wa  is  the 
Bose-Einstein  distribution  for  Ns  photons  with  wavevector  ks  and  frequency  u^, 
and  the  summation  in  Eq.(2.26)  over  the  number  of  photons  can  be  carried  out  ana- 
lytically. We  have  also  used  in  Eq.(2.26);  ((D^ (-t)D^(O))) m = (0^(t)D^(O)))*, 
to  ensure  that  the  transition  rate  Rn  remains  real.  For  the  second  order  processes, 
we  write  the  transition  rates  R22  as 


+00  +00  +00 

R22  = ^EE  / du'  J du  J dt 

VC  V'C'—oo  —00  —00 

<2-27) 

where 


/«, „c  (u',u,t)  = ((*«.  («'.<)' X,<(“.0)))M  (2.28) 

is  the  molecular  polarizability  tensor  TCF  of  the  target,  and 

KcvC  (u>’u’*)  = ((^'C'  (2-29) 
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is  the  electric  field  polarizability  tensor  TCF  of  the  light  source.  Using  Eqs.(2.17- 
2.20)  in  Eq.(2.29),  the  electric  field  polarizability  tensor  TCF  of  the  light  source  can 
also  be  further  decomposed  into  terms  representing  elastic  and  Raman  scattering 
[Appendix  A].  These  terms  when  substituted  back  in  Eq.(2.27),  will  give  transition 
rates  for  second  order  processes  when  thermal  light  is  interacting  with  a large 
molecular  system. 

Equations  (2.22)  and  (2.27)  for  the  transition  rates  Ru  and  R22  are  compu- 
tationally useful,  because  they  do  not  involve  the  summations  over  the  excited 
states  of  the  target,  and  the  corresponding  molecular  TCF  can  be  evaluated  by 
time-dependent  techniques. 

2.3.  Radiation  Field  TCFs  For  Laser  Light 


When  the  light  source  is  a laser,  the  statistical  distribution  of  its  electric  field 
factorizes  into  a distribution  for  the  electric  field  amplitude,  times  a distribution  for 
its  phase  [Louisell,  73;  Sargent,  Scully  and  Lamb,  74;  Loudon,  83].  If  the  laser  is 
operated  below  the  threshold  for  coherent  emission,  the  amplitude  fluctuations  of  a 
given  mode,  which  are  proportional  to  the  mean  number  of  photons  in  that  mode, 
have  an  exponentially  decaying  distribution  about  the  most  probable  amplitude 
(which  occurs  at  zero),  whereas  fluctuations  in  the  phase,  caused  by  spontaneous 
emissions,  are  quite  large  and  therefore  the  phase  can  have  any  value  between  0 
and  27t.  If  the  laser  is  operated  above  the  threshold,  both  the  amplitude  and  the 
phase  have  sharply  peaked  distribution  about  their  most  probable  values,  which  are 
given  by  experimental  conditions.  For  such  a light  source  it  is  necessary  to  include 
in  the  description  the  statistical  distributions  of  both  the  amplitude  and  the  phase  of 
the  electric  field,  and  it  is  advantageous  to  work  in  the  coherent  state  representation 
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[Klauder  and  Sudarshan,  68;  Loudon,  83].  A general  radiation  state  in  the  coherent 
state  representation  is  written  as  | {as}),  with  a ^ | {o3})=  as  | {a,})%,  where 
the  eigenvalue  as  is  a complex  number  whose  square  modulus  gives  the  mean 
number  of  photons  in  the  s mode  of  the  laser.  Using  equations  (2.17-2.20)  for  the 
creation  and  annihilation  terms  of  the  electric  field  operator,  and  using  a single 
mode  coherent  states  | as),  we  write 


The  radiation  field  TCF  of  Eq.(2.22),  in  the  coherent  state  representation,  can  be 
written  as 


Unlike  f^(t)  in  the  photon  number  representation,  cross-terms  in  Eq.(2.31)  are  not 
zero,  because  we  do  not  need  to  conserve  the  number  of  photons  in  the  initial  and 
the  final  radiation  states  | r)  and  |r'),  respectively.  As  we  shall  see  very  shortly, 
however,  it  is  still  possible  to  have  a laser  in  which  the  statistical  distribution  of  the 
phase  is  such  that  it  again  can  make  the  contributions  of  these  cross-terms  equal 
to  zero.  Using  equations  (2.18),  (2.19)  and  (2.30)  in  Eq.(2.31),  we  wnte  in 

a single  mode  coherent  state  representation  as 


(2.30) 


/&  m = ((4_>  «4+)  (°)))„ + ((4° « 4"’  (°>))„ 


(2.31) 


(2.32) 
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where  P(as),  a real  statistical  weight  for  the  complex  amplitude  as,  is  given  by 
experimental  conditions.  The  two  dimensional  integral  f d2a*[...]  extends  over 
the  entire  complex  plane  of  as.  Light  sources  of  experimental  significance  have 
different  statistical  weights  P(as).  In  the  first  and  the  fourth  terms  of  Eq.(2.32),  we 


factors  P(as ) given  in  the  literature  [Schubert  and  Wilhelmi,  86].  Therefore,  the 
radiation  field  TCF  in  Eq.(2.32)  satisfies  the  relation 

and  hence,  the  resulting  transition  rates  of  experimental  significance,  are  real  from 
Eq.(2.32).  As  examples  of  laser  sources,  with  given  statistical  weights,  we  consider 
the  radiation  field  TCFs  for  chaotic  and  coherent  light.  To  compare  Eq.(2.32)  with 
the  equivalant  expression  in  the  photon  number  representation  in  Eq.(2.25),  we 
have  first  considered  the  case  of  chaotic  light,  which  is  the  same  as  that  of  a laser 
operating  below  the  threshold.  In  this  case  the  phase  of  the  light  source  is  randomly 
distributed  between  0 and  27r,  therefore  the  weight  factor  P(a3)  does  not  depend 
on  the  phase  of  a3,  and  hence,  the  contributions  from  the  first  and  the  fourth  terms 
of  the  integrand  in  Eq.(2.32)  are  zero.  Thus,  below  the  threshold  takes  a 

simpler  form 


which  agrees  with  Eq.(2.25)  for  a single  mode,  insofar  as  the  average  is 


have  checked  that  the  statistical  averages  of  a2,  and  a*2  are  real  for  several  weight 


OO 


0 


(2.34) 
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Substituting  Eqs.(2.32)  and  (2.33)  in  Eq.(2.21),  we  can  write  the  first  order  transition 
rate  Rn  for  the  interaction  of  a single  mode  laser  with  a molecular  system,  when 
the  laser  is  operating  below  the  threshold,  as 

OO 

Rcn°“c  = £^7  Ef£  / <“[(*•>  “*("£.0 

( 0 

x2Re(((Dl(t)'D((0)))M)  (2.35) 

(()*££(  0)))m)) 

If  the  laser  is  operating  above  the  threshold,  the  rate  of  transition  involves  instead 
the  statistical  distribution  P(os),  which  depends  on  the  amplitude  | as  |,  as  well 
on  the  phase  = Arg(as),  and  is  given  by 


TtCoherent 

Ku 


Uk. 


+oo 

J dt  J d2as  P(as ) 


^ — OO 

i2 


2he0V 

x [2  |Qa|^  {cos(uj-t)}  - cos(<^<  - 2 <p3) 

x«£{(t),D{(0)»M 


+ e-“.-.'«D£((),D£(  0)))„] 


(2.36) 


This  can  be  rewritten  as 


TyCoherent  nChaotic 

R-ll  ~ KU 


2e0V 


OO 


dt 


x [As  cos (u^t)Re{((D^  (<)f  D{  (0)))m} 


(2.37a) 


+ iBs  sin(u >£t)Im{((D£  (f)f  (0)))ji/}] 
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with 

(2.376) 
(2.37c) 


A s = J d2as  P (as)  |aa|2  cos(2v?s) 
Bs  = J d2as  P(as)|aa|2  sin(2 <pa) 


showing  an  additional  rate  due  to  coherence.  The  term  containing  Ds  can  be 
shown  to  be  zero  for  statistical  distributions  of  experimental  interests  [Schubert 
and  Wilhelmi,  86]. 


CHAPTER  3 

MOLECULAR  TCFs  IN  THE 
ELECTRONIC  ADIABATIC  REPRESENTATION 


In  the  previous  chapter,  transition  rates  for  all  first  and  second  order  pho- 
tointeraction processes  are  expressed  as  integrals  over  time  of  the  product  of  the 
molecular  TCFs  of  the  target  and  the  radiation  field  TCFs  of  the  light  source.  The 
time -dependence  in  the  radiation  field  TCFs  is  expressed  for  many  physically  in- 
teresting light  sources  in  terms  of  exponentials,  oscillating  with  the  frequency  of 
the  light  source  used.  Therefore,  most  of  the  transition  rates  are  finally  written  as 
the  Fourier  transforms  of  the  molecular  TCFs  of  the  target.  In  this  chapter,  we 
describe  general  features  of  these  molecular  TCFs  and  their  Fourier  transforms. 
Indeed,  the  task  of  a general  formalism  should  normally  include  describing  a light 
source  within  a frequency  range  that  extends  from  the  far  infrared  region  on  one 
side,  to  the  ultravoilet  region  on  the  other.  A complete  discussion  of  such  an  ex- 
tended frequency  range,  however,  is  beyond  the  scope  of  this  work,  and  we  will 
describe  mainly  the  interaction  of  visible  or  UV  light  with  large  molecular  systems. 
As  shown  in  table  3-1,  the  light  source  in  this  frequency  range  excites  nuclear 
vibrations,  accompanied  by  electronic  transitions  between  different  states  of  the 
system  [McQuarrie,  83].  Thus,  our  work  involves  mainly  the  electronic  excitations 
of  the  molecular  system,  followed  by  a time-dependent  evolution  according  to  the 
upper  potential  energy  surface,  which  in  some  cases  can  also  cause  the  dissociation 
of  the  system.  We  have  also  considered  nuclear  motions  in  the  ground  electronic 
state,  which  occur  generally  when  the  molecule  is  in  a thermal  bath. 
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Table  3-1 

The  Electromagnetic  Spectrum  And  The  Corresponding  Molecular  Processes 


Region 

Wave  # (/cm) 

Process 

Microwave 

0.033  - 3.300 

Rotations  of  Polyatomic 
Molecules 

Far  Infrared 

3.30  -330 

Rotations  of  small 
molecules 

Infrared 

330  - 3300 

Vibrations  of  flexible 
bonds 

Visible  and  Ultravoilet 

3300  - 330000 

Electronic  transitions 
with  vibrations  of 
flexible  bonds 

[From,  “Quantum  Chemistry”  by  D.  A.  McQuarrie;  83,  pp.438] 


In  section  3.1,  we  introduce  electronic  potential  energy  surfaces  for  the  motion 
of  the  nuclei,  and  write  a general  molecular  TCF  of  the  target  in  the  adiabatic 
representation.  Due  to  the  difficulty  involved  in  computing  a real  time  propagator 
with  realistic  potential  energy  surfaces,  we  describe  a transformation  to  complex 
time  TCFs  in  section  3.2.  Section  3.3  discusses  general  features  of  these  complex 
time  TCFs,  and  writes  them  in  a form  which  is  used  in  the  later  part  of  this  work. 
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3.1.  Electronic  Adiabatic  Representation  for  Molecular  TCFs 


We  consider  the  evaluation  of  a general  positive  time  TCF  for  two  different 
operators  A and  B.  In  case  of  the  first  order  processes,  they  are  the  molecular 
transition  dipole  operators,  and  in  case  of  the  second  order  processes  such  as  elastic 
and  Raman  scattering,  they  become  the  molecular  polarizability  tensor  operators. 
A general  positive  time  molecular  TCF  for  operators  A and  B is  written  as 


hs  (*)  = ((Mt)' 6 (0)))  u (3.1) 

where  ((...))m  is  the  quantal  average  over  the  probability  distribution  of  the  initial 
state  of  the  molecular  system,  and  the  time-dependence  in  the  operator  A is  given  by 


A(t)  = e“^M<i(0)  e+^Mt  (3.2) 


Here  Hm  the  total  Hamiltonian  for  the  extended  molecular  system  involves  the 
nuclear  as  well  as  the  electronic  motions  in  the  system,  and  is  written  as 


Hm  = - £ 


3N-6  *2 


d2  h2  d 2 . 

-E—  — + ^(r’R) 


^ 2MndR2  ^2 mdr2 

n=l  «=1  1 


(3.3) 


with  r = (ri,r2,r3, ...)  is  a set  of  s generalized  coordinates  for  the  motion  of 
electrons,  and  R = (Ri,  R2,  R3,  •••)  is  a set  of  3 N -6  generalized  coordinates  for  the 
vibrational  motions  of  the  nuclei.  We  use  the  Bom-Oppenheimer  approximation  to 
separate  the  electronic  and  nuclear  motions  in  the  target.  The  description  is  suitable 
for  the  ground  and  first  few  electronically  excited  states  of  the  target,  where  such 
a separation  in  the  electronic  and  nuclear  motion  generally  yields  reasonably  good 
results.  In  the  Bom-Oppenheimer  approximation  [Balint-Kurti,  75]  the  electronic 
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wave  functions  are  the  solutions  of  the  following  equations 

I"  E J?  + V<r'  R'l  ^ R>  = Ei(R)  R)  (3-4) 

1 = 1 1 

This  equation  is  parametric  in  the  nuclear  variables  R,  and  {Ej(R)}  is  a set  of 
electronic  potential  energy  surfaces  for  the  motion  of  the  nuclei.  The  total  initial 
wave  function  of  the  target  is  then  written  as 

I*/, (R)>  = Vi(R)  |*i(R)>  (3.5) 

where  |<^(R)),  the  initial  electronic  state,  depends  parametrically  on  the  nuclear 
variables  R,  |i/>/)  is  the  nuclear  vibrational  wave  function  on  the  ith  electronic  po- 
tential energy  surface,  and  the  explicit  dependence  of  the  electronic  wave  functions 
|<^j(R))  over  the  electronic  coordinates  r is  suppressed,  because  they  are  integrated 
out  in  the  next  step.  Then,  if  Wj  is  the  probability  of  finding  the  nuclear  vibra- 
tional wave  function  in  the  vibrational  state  of  the  ith  electronic  potential 
energy  surface,  the  TCF  /^(f)  as  defined  in  Eq.(3.1)  can  be  written  as 

h a(0  = E w,  (<S,(R)  | (*  | (A) 

il  (3.6) 

e-i*»‘(B)|,h>|*(R)) 

Introducing  a complete  basis  set  in  the  coordinate  representation 

E |R)(R|V»>I«R))Wj(R)I(0j|R)(R|  = 1 . (3.7) 

jJ 

and  using  the  adiabatic  approximation  for  electronic  potential  energy  surfaces 
[Olson  and  Micha,  82]  expressed  as 

£j#l*(R)«  |<MR)>Wi(R,PR)  (3.8) 

we  write  Eq.(3.6)  as 

fAj(t)  = E E W‘  f dR  (A(R))ij 

ij  I J 


(3.9) 
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where  ?fi(R,PR)  is  the  Hamiltonian  for  the  nuclear  motion  on  the  ith  potential 
energy  surface,  and  Pr  is  the  set  of  the  momentum  operators  conjugate  to  the 
operators  in  the  set  R.  In  the  time -dependent  formulation,  we  interpret  the  TCF 
f a ^(0  *n  Eq.(3.9)  as  follows:  at  an  initial  time,  0,  the  amplitude  of  the 
nuclear  motion  on  the  jth  electronic  surface  is  given  by  (-B)j,  | t pj)  . At  a 
later  time,  t1  > 0,  this  amplitude  propagates  on  the  excited  electronic  surface 
with  the  nuclear  Hamiltonian  Tij  = 7fj(R,PR).  Finally,  at  the  required  time 
t'  = t,  the  time  evolved  amplitude  on  the  jth  electronic  surface  is  overlapped 
with  the  complex  conjugate  of  | tpi).  In  case  the  system  has 

initially  been  prepared  in  a pure  vibrational  state  | ?/»/)  of  the  ith  potential  energy 
surface,  the  initial  probability  distribution  function  Wj  becomes  unity,  and  the 
TCF  in  Eq.(3.9)  describes  a general  light-molecule  interaction  process  involving 
two  different  electronic  states  i and  j.  If  the  system,  however,  is  not  initially  in  a 
pure  state  | ipj),  then  the  finite  temperature  formalism,  in  which  the  initial  state  of 
the  system  has  been  prepared  in  a thermal  equilibrium,  can  be  introduced  simply 
by  selecting  a suitable  probability  distribution  function  Wj.  Such  a probability 
distribution  function  is  usually  described  by  the  Boltzmann  distribution  expressed 
as  Wj  = e(  ^r— , where  /?  is  the  inverse  temperature,  ksT,  at  which  the  initial 
state  of  the  system  has  been  prepared,  and  £j,  the  stationary  eigen  state  on  the  ith 
potential  energy  surface,  is  given  by 

Wi)  = w i)  (3.10) 

For  a two  surface  problem  using  Eq.(3.10)  in  Eq.(3.9),  the  TCF  f ^ ^(<)  is  expressed 
in  the  following  equivalent  forms: 
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dR  f dR'  (R|  e+5^(t+i/?h)  |R')  A(R') 


(R'|e-5^t|R)B(R)  (3.116) 


= ^wIe+i£'t{ti\(A)ije-Mt(B)ji : |V>/>  (3.11c) 


/ 


where  ?fj  = ftj(R,  Pr)  is  the  nuclear  Hamiltonian  operator  on  the  jth  potential  en- 
ergy surface,  and  the  !>[...]  in  Eq.(3.11a)  is  to  be  evaluated  over  the  initial  nuclear 
vibrational  state  \xpj)  of  the  system.  Equation  (3.11b)  rewrites  this  trace  and  all  the 
operators  in  the  coordinate  representation,  which  is  useful  especially  if  the  initial 


and  the  j*  electronic  surfaces  are  not  difficult  to  obtain.  The  form  in  Eq.(3.1  lc)  is 
used  when  the  initial  nuclear  vibrational  state  | xpj)  and  the  corresponding  station- 
ary state  eigenenergy  E\  are  not  difficult  to  obtain  by  time-independent  techniques, 
while  the  important  effects,  such  as  dissociation  and  rearrangement,  are  occuring 
due  to  the  time-dependent  dynamics  on  the  excited  surface. 

Equations  (3.11a),  (3.11b),  and  (3.11c)  present  a general  two-operator  TCF 
in  the  adiabatic  representation,  where  (A)tJ  and  (B)ji  are  the  transition  moment 
operators  between  two  different  electronic  states  i and  j.  If  the  frequency  range 
of  the  light  source  is  in  the  infrared  region,  then  it  excites  nuclear  vibrations  and 
rotations  only  on  the  ground  electronic  surface  of  the  system.  In  our  work,  this 
corresponds  to  the  case  j = i,  with  i refering  to  the  ground  electronic  state  of  the 
system.  Transition  moments  ( A)tt  and  (B)a  in  this  case  become  the  net  dipole, 
or  the  net  polarization  of  the  system  in  its  ground  electronic  state,  and  the  TCF  as 


nuclear  vibrational  state  of  the  system  is  not  known,  and  the  propagators  on  the 
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expressed  in  Eq.(3.11a)  acquires  a symmetric  form  which  can  also  be  evaluated  in 
our  time -dependent  formulation. 


3.2.  Complex  Time  TCFs  and  their  Fourier  Transforms 


The  usual  route  for  computing  the  transition  rates  for  various  different  light- 
molecule  interaction  processes  that  emerges  from  the  discussion  till  now,  is  to  first 
compute  the  time  dependence  in  the  TCF  as  defined  in  Eq.(3.11),  and  then  take 
its  Fourier  transform  to  get  the  required  transition  rate.  The  route  as  suggested, 
however,  is  not  very  easy  to  follow  for  any  physically  interesting  problem  involving 
realistic  potential  energy  surfaces.  This  is  primarily  due  to  the  presence  of  a real 
time  excited  potential  propagator  in  the  expression  for  the  TCF  in  Eqs.(3.11a, 
b,  c).  The  form  of  this  propagator  is  more  apparent  in  Eq.(3.11b),  where  it  is 
expressed  in  the  coordinate  representation.  Computations  of  the  matrix  elements 
of  a pure  real  time  propagator,  involving  realistic  potentials,  is  a very  difficult 
problem  in  itself,  so  that  to  our  knowledge,  till  now,  there  is  no  known  numerical 
or  analytical  method  that  even  tries  to  address  this  problem  for  a completely  general 
potential.  This  is  also  the  main  reason  that  path-integral  methods  [Feynmann 
and  Hibbs,  65;  Schulman,  81]  for  evaluating  a propagator  in  real  time  dynamics, 
though  very  simple  to  understand  and  physically  appealing  to  the  intuition,  have 
remained  basically  a formal  device,  and  have  not  been  used  in  problems  involving 
realistic  potentials.  Recendy,  there  has  been  a resurgence  of  interests  in  these 
methods  [Miller,  Schwartz  and  Tromp,  83;  Behrman,  Jongeward  and  Wolynes,  83; 
Thirumalai  and  Beme,  83;  Doll,  84],  and  significent  advances  are  made  towards 
developing  numerical  techniques  for  the  evaluations  of  rapidly  oscillating  integrals 
in  many  dimensions  [Chang  and  Miller,  87;  Doll,  Coalson  and  Freeman,  87]. 
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Alternatively,  there  have  been  developments  such  as  time-independent  quantum 
mechanical  coupled  channel  methods  [Shapiro  and  Bersohn,  80;  Heather  and 
Light,  83],  the  time-dependent  quantum  mechanical  wave-packet  methods  [Heller, 
81;  Mowrey  and  Kouri,  85;  Sawada  and  Metiu,  86],  and  a host  of  other  semi- 
classical  methods  [Marcus,  72;  Miller,  74;  Mukamel  and  Jortner;  76,  Micha  83; 
Billing  and  Jolicard,  84;  Gray  and  Child,  84;  Henriksen,  85]  which  are  generally 
based  on  running  classical  trajectories  for  the  nuclear  motions,  while  the  electronic 
transitions  are  treated  quantum  mechanically.  Each  of  these  methods  have  their  own 
advantages  and  disadvantages,  and  have  been  applied  to  the  electronic  absorption 
and  photodissociation  problems  with  varying  degrees  of  success.  Finite  element  grid 
methods  which  are  based  on  the  path-integral  approach  to  computing  a propagator 
in  the  coordinate  or  the  phase  spaces  [Thirumalai,  Bruskin  and  Berne,  83;  Coalson, 
85;  Jackson  and  Metiu,  85],  and  are  not  restricted  to  any  specific  form  for  the 
potential,  have  not  been  applied  to  any  of  these  light-molecule  interaction  processes 
involving  realistic  potential  energy  surfaces  in  more  than  one  dimension.  This  can 
be  explained  very  briefly  by  considering  the  coordinate  space  matrix  elements  of 
a real  time  propagator  after  a time  interval  2e  as 


where  H is  a general  one  dimensional  time-independent  Hamiltonian  along  the 
variable  A^.For  e — *•  0,  a typical  expression  for  a short  time  approximation  for  the 
matrix  elements  appearing  in  the  integrand  of  Eq.(3.12)  [Schulman,  81],  is  written 
as 


+ 00 


(3.12) 


(X"\e~tH(\X) 
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(313> 

£-00  -.|(y(X)  + V(X'))} 

The  first  term  in  the  exponent,  which  corresponds  to  the  short  time  approximation 
for  the  kinetic  energy  term  in  the  propagator,  increases  as  the  square  of  the 
difference  (A  — A').  The  oscillating  nature  of  the  exponentials  therefore  increases 
rapidly  as  one  considers  larger  and  larger  values  of  (A  —X'),  whereas  the  magnitude 
of  these  rapidly  increasing  oscillations  remain  constant  over  the  entire  coordinate 
space.  All  of  the  numerical  techniques  for  computing  the  integrals  in  Eq.(3.12), 
which  are  based  on  approximating  the  limits  of  the  infinite  integrals  in  Eq.(3.12) 
to  be  within  some  large  finite  region  of  the  coordinate  space,  therefore  do  not 
converge  to  any  physically  meaningful  result.  One  approach  to  avoid  this  problem 
is  to  use  a transformation  which  shifts  the  problem  of  computing  the  real  time  TCF 
to  another  problem,  in  which  instead  one  computes  the  same  TCF  for  complex 
times,  and  in  the  end,  after  computing  the  Fourier  transform  of  this  complex  time 
TCF,  one  uses  some  suitably  defined  inverse  transformation  to  extract  the  required 
physical  quantity  of  interest.  The  real  to  the  complex  time  transformation  of  the 
computations  in  Eq.  (3.12),  the  real  short  time  e — > oo  is  replaced  by  an  equivalent 
complex  short  time  k — ► oo,  where  k has  a fixed  non-zero  negative  imaginary  part. 
As  a result,  the  rapidly  oscillating  exponentials  in  Eq.(3.12),  which  have  a uniform 
maginutde  in  the  entire  coordinate  space,  are  replaced  by  the  rapidly  oscillating 
exponentials  the  magnitude  of  which  also  decay  rapidly,  as  one  considers  larger 
and  larger  values  (A  — A').  Therefore,  the  numerical  evaluation  of  the  infinite 
integrals,  within  sufficiently  large  finite  limits  in  the  coordinate  space  in  Eq.(3.12), 
does  not  cause  any  problem. 
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Next,  we  describe  the  transformation  from  real  to  the  complex  time  for  a 
molecular  TCF  as  written  in  Eq.(3.11),  and  derive  the  required  inverse 

transformation  that  gives  us  the  required  Fourier  transform  of  the  original  TCF.  If 
initially  the  molecular  system  is  in  a pure  nuclear  vibrational  state  |V>/)  of  the  ith 
potential  energy  surface,  f ^ t ) of  Eq.(3.11c)  can  then  be  written  as 

fA6  (<)  = (3.14) 

I,J 

where  | tpj)  is  the  final  nuclear  vibrational  state  on  the  jth  potential  energy  surface, 
and  ujj  = (£]  — £i)/h.  The  Fourier  transform  of  Eq.(3.14),  defined  as 

+oo 

= W (313) 

— oo 

gives  the  familiar  Golden  rule  result 

- “»)  (3.16) 

I yJ 

where  the  Dirac  6 function  conserves  the  energy,  by  requiring  that  the  energy  hujj 
of  the  electronic  transition  between  the  initial  and  the  final  nuclear  vibrational  states 
| rpj)  and  |t/>j)  is  equal  to  the  energy  hu  given  by  the  light  source.  The  Fourier 
transform  as  defined  in  equations  (3.15),  requires  that  the  computations 

are  done  along  the  real  time  axis  in  the  complex  time  plane.  The  transform  to  an 
axis  in  the  complex  time  plane,  which  is  parallel  to  the  real  time  axis,  but  with  a 
non-zero  constant  imaginary  part,  is  achieved  simply  by  replacing  in  Eq.(3.14)  the 
real  time  t by  a complex  time  r = t — ihu,  i.e  we  use  = /^(t) 
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hi>w  = <l>j)  (3.17) 

(VvK  £).,.#/) 

where  u is  a parameter  which  always  remains  real.  Using  the  definition  in  Eq.(3.15) 
for  the  Fourier  transform,  we  can  see  that  the  Fourier  transform  of  the  complex 
time  TCF  f ^ g(t)  is  given  by 


hi  (")  = (3.18) 

* u 

( ipj\(B)ji\ipi ) 6(u-uji) 

From  equations  (3.16),  (3.18),  and  using  the  property  of  the  Dirac  6 function 
requiring  ujjj  = u,  it  follows  that  which  defines  the  transition  rates  for 

different  light-molecule  interaction  processes  in  the  previous  chapter,  is  given  by 

(«’)  = '"**/«(")  (3-19) 

where  is  the  Fourier  transform  of  the  complex  time  TCF  as  defined 

in  Eq.(3.17).  For  a constant  u this  transform  is  taken  along  an  axis,  which  is  parallel 
to  the  real  time  axis  but  has  a non-zero  constant  imaginary  part.  In  case  the  initial 
state  of  the  molecular  system  is  not  a pure  vibrational  state,  but  is  at  thermal 
equilibrium  at  the  inverse  temperature  we  use  Eq.(3.11b)  of  the  previous 
section  and  steps  similar  to  Eqs.(3.14-3.19),  to  show  that  if  the  complex  time  TCF 
for  a finite  temperature  two-surface  photoexcitation  problem  is  given  by 
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/,4M=^Tr 


e+*(t+,M)w.  / 

4)  e-tt*-™)*'  L 

b)  | 

V 

) ij  V 

/>«J 

(3.20) 


with  ft  = ft  — u,  and  u as  the  artificially  introduced  parameter,  then  f ^ g(u>)= 
euhu>f^  where  /^(w)  is  the  Fourier  transform  of  Eq.(3.20)  as  defined  in 
Eq.(3.15).  Some  of  the  other  works  along  these  lines  treat  the  artificially  introduced 
parameter  u as  temperature;  the  computations  of  a complex  time  propagators  are 
then  suitably  referred  to  as  real  time  finite  temperature  computations.  In  case  of  an 
actual  finite  temperature  computation,  due  to  a thermal  distribution  for  the  initial 
state  of  the  system,  the  requirement  of  a positive  ft  in  the  argument  of  the  lower 
surface  propagator  in  Eq.(3.20)  imposes  therefore  an  upper  limit  on  the  choice  of 
the  parameter  u.  In  our  work,  however,  we  do  not  use  time-dependent  computations 
for  the  propagator  on  the  lower  potential  surface;  therefore,  we  have  not  needed  any 
constraint  on  the  size  of  the  parameter  u imposed  by  the  size  of  the  natural  ft  for  the 
problem.  We  continue  to  call  u simply  an  artificially  introduced  parameter,  the  size 
of  which  depends  only  upon  the  following  competing  factors.  Increasing  the  value 
of  u increases  the  damping  of  the  upper  surface  oscillations  and  therefore  increases 
the  numerical  stability  of  the  dynamical  computations.  Increasing  the  value  of  u, 
on  the  other  hand,  also  decreases  the  size  of  the  numbers  we  are  dealing  with  in 
the  computations;  therefore,  beyond  a certain  value,  they  start  to  affect  the  overall 
accuracy  of  the  computations.  For  each  individual  problem,  one  has  to  choose  a 
compromise  between  these  two  competing  factors  and  select  a suitable  range,  such 
that  u is  sufficiently  large  for  an  effective  damping  for  the  upper  surface  oscillations, 
and  at  the  same  time  sufficiently  small  so  that  the  numbers  in  the  computations  are 
not  too  small  to  affect  the  overall  accuracy.  In  the  end,  the  accuracy  of  a particular 
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computation  is  measured  by  checking  the  stability  of  a computation  with  respect 
to  the  variations  in  the  value  of  u. 

The  transformation  of  a computation  along  the  real  time  axis  to  another  com- 
putation along  a parallel  axis,  but  with  a non-zero  imaginary  component  u,  is  not 
unique.  There  are  other  transformations,  such  as  rotating  the  computational  axis 
in  the  complex  time  plane,  that  can  achieve  similar  results.  In  fact,  some  of  them 
may  later  prove  to  be  more  efficient  than  the  transformation  used  in  this  work.  We 
discuss  this  further  in  the  later  chapters. 

Lastly,  the  property  of  the  <5  function  for  the  conservation  of  energy  in 
Eqs.(3.16)  and  (3.18),  which  has  been  used  in  Eqs.(3.17)  and  (3.20)  to  get  the 
inverse  transformation  between  f ^ ^(u)  and  f ^ g(u),  holds  true  only  for  a sharp 
spectra.  Other  cases,  where  one  has  to  incorporate  also  a phenomenological  line 
shape  in  the  problem,  such  as  molecular  line-broadening  effects,  we  can  use  mod- 
ifications of  the  similar  transformation  [Coalson,  85]  to  achieve  similar  results. 

3.3.  Properties  of  the  Complex  Time  TCFs  and  their  Fourier  transforms. 

Some  properties  of  the  complex  time  TCF  f ^ g(t),  which  we  use  later  in  this 
work,  follow  from  its  explicit  expression  which  we  write  once  again  as 


/.«(<)  = £ w,U, 


e+**'r  L 

i)  e-Wr(l 

9) 

V 

) ij  V 

*i'i 


(3.21) 


The  initial  and  the  excited  surface  Hamiltonians  7i[  and  Hj,  respectively,  are  the 
operators  over  the  nuclear  variables  {R,Pr}  of  an  extended  molecular  system. 
Sometimes  these  Hamiltonians  may  contain  constant  terms  like  hA,  and  fiAj 
respectively,  such  that  the  difference  hAXJ-  h(Aj  - A,)  may  correspond  simply  to 
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a vertical  displacement  of  the  upper  surface  with  respect  to  the  lower  one.  These 
constant  terms  affect  only  the  position  of  the  line  shape  function  in  the  frequency 
space.  Therefore,  they  can  be  taken  care  of,  very  easily,  in  the  definition  of  the 
Fourier  transform  f ^ g(u>).  If  the  upper  and  the  lower  surface  Hamiltonians  are 
written  as  Hj=  H}  - hAj  and  7Y,=  Ht  - TiAt,  respectively,  f ^ g in  Eq.(3.22) 
can  then  be  written  as 

W = (0  <3-22) 

with  its  Fourier  transform  given  by 

= • <3-23) 
The  shifted  Fourier  transform  F^g(u  — A}t)  is  given  by 

+oo 

= j dt ‘“'hi W <3-24) 

— oo 

with  J = u — A ji,  and  the  shifted  TCF  F ^ g(t)  is  written  in  a compact  form  as 


I N 


(A)  e~'^T  (l 

b) 

V ) ij  V 
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(3.25) 


Here  we  have  used  H,\ipi)  = Ej\ipi),  and  the  frequency  operator  Clj(Ej)  which 
describes  the  time-evolution  is  given  by 


(g>  - Ei) 


h 


(3.26) 
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where  Hj  is  an  operator  over  only  the  nuclear  variables  {R,Pr}  for  the  jth 
electronic  surface.  For  two  different  operators  A and  B,  we  define  a symmetric 
form  for  the  TCF  as 


= + <3-26> 

which  is  invariant  under  the  interchange  of  operators  A and  B.  Replacing  t by 
—t  in  Eq.(3.25) 


= £**'>  (*  I (^)jie+,A'(R)r‘  (*),, 


xPl 


(3.28) 


we  see  that  F ^ g(-t)  = F ^ f )*.  Now  considering  the  Fourier  transform  F^  ^(u>) 

of  Eq.(3.25)  as 


( 0 
A 

+oo  \ 

J dt'e"<'FAt(t)  + 

/ FAb  (0 

~oo 

o J 

replacing  in  the  first  term  t'  by  -t',  and  using  the  property  F ^ g(-t')  = Fg 
we  get 

+oo 

hi  M = ^ / *'  (')•  + (<)}  (3.30) 

0 

For  B = A this  gives 

+oo 

FaaM  = 1 J dt'Re{e^‘'FAA(t)} 

0 


(3.31) 
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and  for  B ^ A,  we  use  the  symmetrized  form,  which  is  invariant  to  the  interchange 
of  the  operators  A and  B,  to  write 


+oo 

= ; / dt'Re{c^F-A6(t)}  (3.31) 

0 

Equations  (3.31)  and  (3.32)  for  the  Fourier  transforms  F^g( w)  and  Fs^^(u>), 
respectively,  ensure  that  the  physical  properties  represented  by  these  Fourier  trans- 
forms remain  real.  Furthermore,  for  B ^ A,  we  can  show  that  the  symmetrized 
TCF  Fs-  -(f)can  also  be  written  as 


fAB  M = \ {%)  (*)  - PXX  (*)  - PBB  (<)}  (3  33) 

where  F^  ^(t)  and  F^  ^(t)  are  the  autocorrelation  functions  correspond- 

ing to  the  operators  C,  A and  B,  respectively,  and  C = A + B.  Therefore,  for  all 
physically  interesting  light-molecule  interaction  processes,  from  now  on,  we  con- 
sider the  evaluation  and  uses  of  the  complex  time  autocorrelation  function  F^^(t), 
which  we  again  write  in  a compact  form  as 


PAA  « = E W>  (Ti.  (°)  lT'.  (’•))  (3.34) 

I 

where  the  time-dependent  amplitude  | T^,(f)),  which  is  evolving  according  to 
the  frequency  operator  is  given  by  |Ty,(f))  = exp(— tftjf)|Tj,-(0))  with 

|Tj-(0)  = (A)ji\ipj).  The  overlap  (Tj,(0)  | Tj-(r)),  completely  describe  the  time- 
dependence  in  the  TCF  F ^ ^(f)  for  a complex  time  r = t — iuh.  To  compute  a 
vibrational  infrared  or  a rotational  microwave  spectra,  we  merely  need  to  replace 
the  frequency  operator  Qj  by  a different  frequency  operator  Q,=  (Ht  — Ei)/fi , 
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and  redefine  the  transition  dipole  moment  operator  (A)ji  as  (.4),,;  the  total  dipole 
operator  of  the  system  in  its  ith  electronic  state. 

The  physical  interpretation  of  Eq.(3.34)  is  as  follows:  at  an  initial  real  time 
t = 0 and  the  complex  time  r = — iuh,  the  amplitude  | Tj,-(— iuh))  is  almost 
identical  in  shape,  but  of  smaller  magnitude  in  comparison  to  | Yj-(O)).  For  later 
times  t > 0,  the  propagation  of  the  amplitude  occurs  according  to  the  frequency 
operator  Clj,  and  a fall  off  in  the  value  of  the  overlap  (Tj,(0)  | T j,(r))  is  observed, 
as  the  amplitude  | T j,(r))  moves  away  from  its  initial  position  and  begins  to  change 
its  shape.  The  short  time  behavior  of  the  overlap  primarily  decides  the  position  and 
the  shape  of  the  overall  line  shape  in  the  high  frequency  region,  whereas  the  long 
time  behavior  of  the  overlap  decides  its  shape  in  the  low  frequency  region.  If  the 
overlap  decays  to  zero  at  intermediate  times  and  shows  some  recurrences  in  the  long 
time  behavior,  one  finds  fine  structures  overlaying  the  otherwise  smooth  spectra. 


CHAPTER  4 

FACTORIZATION  OF  TCFs  INTO 
PRIMARY  AND  SECONDARY  REGION  TCFs 


We  consider  the  evaluation  of  a general  complex  time  dipole-dipole  TCF 
Ffo  p(t)  for  light-molecule  interaction  processes,  involving  electronic  transitions 
between  two  potential  energy  surfaces.  In  particular,  we  are  interested  in  application 
to  the  photoexcitation  and  dissociation  dynamics  of  extended  molecular  systems. 
The  usual  problem  one  encounters  in  dealing  with  extended  molecular  systems 
is  the  large  number  of  degrees  of  freedom  involved.  If  we  were  to  apply  any 
quantum  molecular  dynamics  method  for  computing  a molecular  dipole  TCF  of  an 
extended  system  with  realistic  potential  energy  surfaces,  the  first  and  the  foremost 
consideration  we  would  have  to  think  of  is  the  computational  efficacy  of  the  method 
for  a multivariable  problem.  It  has  been  observed  that  the  computational  efforts 
for  most  of  the  methods,  where  there  is  more  than  one  degree  of  freedom  in  the 
problem,  increase  as  the  square  of  the  number  of  degrees  of  freedom  involved. 
Therefore,  all  the  methods  that  try  to  compute  the  dynamics,  simultaneously  along 
all  the  coupled  degrees  of  freedom,  become  more  and  more  impractical  as  the 
number  of  degrees  of  freedom  start  to  increase  from  more  than  one  or  two.  To 
avoid  this  problem,  we  treat  the  dynamics  of  an  extended  molecular  system  in 
the  spirit  of  the  well-known  mean  field  approximation,  which  factorizes  the  full 
dynamics  of  an  extended  molecular  system  into  many  single  variable  problems, 
with  time-dependent  self-consistent  couplings  between  the  motions  along  these 
uncoupled  variables.  As  a result  of  such  factorization  the  computational  efforts 
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in  a multivariable  dynamical  problem  do  not  scale  as  the  square  of  the  number  of 
the  variables  involved. 

In  section  4.1,  we  describe  the  sufficient  conditions  for,  a)  defining  fast  and 
localized  photoexcitation  processes  in  a large  molecular  system,  interacting  with 
visible  or  UV  light,  and  b)  factorizing  the  molecular  dipole  TCF  for  these  cases  into 
self-consistent  primary  and  secondary  region  TCFs.  Section  4.2  describes  a time- 
dependent  variational  principle  (TDVP)  for  the  amplitudes  | T7(<)),  introduced 
in  the  previous  chapter.  Section  4.3  derives  a norm  conserving  time-dependent 
self-consistent  field  (TDSCF)  approximation  for  factoring  these  amplitudes  into 
self-consistent  primary  and  secondary  region  factors,  and  finally  in  section  4.4,  we 
use  the  results  of  sections  4.1,  4.2,  and  4.3  to  factor  a dipole-dipole  TCF 
into  self-consistent  primary  and  secondary  region  TCFs. 

4.1.  Sufficient  Conditions  for  the  Factorization  of  the  Amplitudes 


In  general,  light-molecule  interaction  processes  in  extended  molecular  systems 
are  classified  into  two  broad  groups.  First,  the  intramolecular  process,  in  which 
a single  molecule  involving  several  vibrational  and  rotational  degrees  of  freedom 
absorbs  light  in  a given  frequency  range,  depending  upon  which  it  dynamically 
evolves  on  the  ground,  or  excited  electronic  surface.  Second,  the  intermolecular 
process,  in  which  an  extended  molecular  system  absorbs  light  in  a given  frequency 
range,  followed  by  a dynamical  evolution  of  the  system  in  a given  electronic  state. 
For  both  cases,  if  the  absorption  of  light  in  the  extended  molecular  system  is  a fast 
process,  we  can  define  a sufficient  condition  to  describe  a localized  region  in  the 
large  molecular  system,  such  that  this  region  plays  a more  dominant  role  in  the 
light-molecule  interaction  as  compared  to  the  rest  of  the  system.  If  the  frequency 
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of  the  light  source  is  in  the  visible  or  UV  region,  the  sufficient  condition  that 
defines  such  a localized  region,  is  that  the  transition  dipole  operator  (D( R))y  of 
the  system,  is  an  operator  over  the  nuclear  variables  of  the  localized  region  only. 
In  other  words,  the  transition  dipole  operator,  for  a localized  electronic  excitation 
process  in  a large  molecular  system,  is  written  as 

(D(  R))ij  = (D(X))ij  (4-1) 

where  R = {X,  Q}  is  a set  of  generalized  nuclear  variables  for  the  entire  molecular 
system,  with  X = {Xi,X2,X3,...}  a subset  that  defines  the  variables  of  the 
the  localized  region  within  which  the  important  dynamics  takes  place,  and  Q — 
{Qi,Q2,Q3,...}  the  subset  that  contains  all  the  remaining  variables  of  the  system. 
We  use  “primary”  region  for  the  region  covered  by  the  elements  of  the  subset 
X,  and  “secondary”  region  for  the  region  covered  by  the  elements  of  the  subset 
Q.  With  the  definitions  of  the  primary  and  the  secondary  regions  of  an  extended 
molecular  system  at  an  initial  time  t = 0,  in  place,  the  sufficient  condition  for  the 
factorization  of  the  amplitude  | T7(0))  into  factors  corresponding  to  the  primary 
and  secondary  region  dynamics  for  a later  time  t > 0,  is  described  by  writing  the 
initial  nuclear  vibrational  state  | ipj)  on  the  ith  potential  surface  as 

i =i  r,)  i (4.2) 

where  p corresponds  to  the  primary  region  dynamics  and  s to  the  secondary  region 
dynamics.  The  stationary  state  eigenvalue  Ej,  corresponding  to  the  sufficient 
condition  in  Eq.(4.2),  is  then  written  as  Ej  = E’j+Ej  + EVjS , where  the  last 
term  is  due  to  the  coupling  between  the  primary  and  the  secondary  regions  on  the 
ith  potential  energy  surface.  Furthermore,  it  is  not  difficult  to  see  that  the  sufficient 
condition  in  Eq.(4.2),  for  the  factorization  of  the  amplitude  into  self-consistent 
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primary  and  secondary  region  factors,  is  the  validity  of  the  time-independent  Hartree 
(TH)  approximation  for  the  stationary  state  vibrational  wavefunction  on  the  ith 
potential  surface  [Bowman,  86].  Given  the  sufficient  condition  in  Eq.(4.1)  for  the 
definition  of  a localized  region  in  an  extended  molecular  system  is  satisfied,  it  is 
not  difficult  to  find  a time-independent  Hartree  approximation  to  the  stationary  state 
eigenfunction  | ipj)  on  the  ith  potential  energy  surface,  such  that 

|T7(0)}=  |f'(0))h'(0))  (4.3) 

where  | £7(0))  = (D(X)) y | is  the  initial  value  of  the  amplitude  for  the 
primary  region  dynamics,  and  | Tjj( 0))  =|  rpj)  is  the  initial  value  of  the  amplitude 
for  the  secondary  region  dynamics.  As  we  shall  see  later,  equations  (4.1)  and 
(4.2)  completely  define  the  sufficient  conditions  for  the  factorization  of  the  time- 
dependent  amplitude  into  self-consistent  motions  corresponding  to  the  primary  and 
secondary  region  dynamics.  Therefore,  these  are  also  the  sufficient  conditions 
for  factoring  the  molecular  dipole  TCF  for  an  extended  system  into  primary  and 
secondary  region  TCF  factors. 

If  the  frequency  of  the  light  source  is  in  the  infrared  region,  the  sufficient 
condition  for  describing  a localized  process  in  an  extended  system  in  Eq.(4.1)  can 
be  generalized  by  using  a series  expansion  for  the  total  dipole  moment  of  the  system 
in  a single  electronic  state  to  write  it  as 

(i>(R))i,  = £(£>' (XJJijPKQ)  (4.4) 

1 

where  (D\X)u  is  the  net  dipole  moment  of  the  primary  region  when  the  molecular 
system  in  the  ith  electronic  state,  and  P/(Q)  is  a polynomial  in  the  variables  of  the 
secondary  region.  The  description  in  the  later  chapters  of  this  work,  with  slight 
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modification,  can  also  incorporate  this  generalization  of  the  sufficient  conditions  in 
Eqs.(4.1)  and  Eq.(4.2),  to  include  the  interaction  of  infrared  light  with  extended 
molecular  systems. 


4.2.  Norm  Conserving  Time-dependent  Variational  Principle  (TDVP) 


The  time-dependent  amplitude  |T(f))  which  defines  the  time-dependence  of 
the  total  TCF  in  Eq.(3.34)  is  a solution  of  the  following  equation  of  motion 

(;3,-fi)|T(f))  = 0 (4.5) 

where  the  initial  condition  |To)  = |T(0)),  defines  the  normalization  TV  at  an  initial 
time  t = 0,  as  N = (To  | To).  For  a two  surface  problem,  we  have  dropped  in 
Eq.(4.5)  specific  references  to  the  subscript  j in  the  frequency  operator  Clj  and  to 
the  subscript  I in  the  amplitude  | Y7).  When  the  exact  solution  |T)  = |T(<)) 
of  equation  (4.5)  is  complicated,  it  is  sometime  desirable  to  replace  it  by  an 
approximate  solution  |T)  = |Y(<))  with  the  initial  conditions  |To)  = |T(0))  at 
t — 0.  We  construct  an  approximate  solution  to  Eq.(4.5)  by  using  a time-dependent 
variational  principle  (TDVP),  the  general  form  of  which  is  similar  to  those  given 
by  Dirac  [Dirac,  26,  30]  and  Frenkel  [Frenkel,  34]  in  the  earlier  days,  and  by 
Mclachlan  [McLachlan,  64]  in  later  modifications,  see  for  example  [Langhhoff, 
Epstein  and  Karplus,  72]  for  a review.  In  our  work,  however,  we  also  impose  the 
conservation  of  the  total  norm  N by  using  the  well  known  Lagrange  undetermined 
multiplier  technique. 

The  variational  functional  is  constructed  in  steps,  by  first  introducing  the 


Lagrangian 


44 


L [T,  Y*]  = ^ ((T|PT)  + (£t|t)) 


(4.6) 


where  V = idt  — &,  and  the  form  of  the  Lagrangian  in  Eq.(4.6)  is  invariant  under 
time  reversal.  From  this,  the  action  between  the  initial  and  the  final  times  follows  as 


where  in  our  notation  the  subscript  A in  the  approximate  solution  |T*)  refers  to  the 
undetermined  Lagrange  multiplier  A,  which  is  solved  by  using  the  constraint  that 
the  norm  N defined  as  N = (T \(t)  | T \(t))  is  a conserved  quantity.  We  solve 
equation  6S'[T^,  T>]  = 0 with  fixed  boundary  conditions,  and  treat  the  variations 
and  6Y\  independendy  of  each  other.  We  choose  the  variation  SY\  = 0 and 
allow  for  a to  get 


(4.7) 


and  finally  the  norm  conserving  functional  for  Eq.(4.5)  is  written  as 


(4.8) 


(«-a|PTj  + ATa)  = 0 


(4.9) 


which  is  equivalent  to 


VTX  + ATa  = 0 


(4.10) 


The  subscript  A in  the  approximate  solution  Y\  refers  to  the  additional  constraint, 
which  needs  to  be  satisfied  solving  for  A.  Using  a transformation  to  remove  A from 
Eq.(4.10),  we  write  the  trial  solution  as 


45 


where  $a  satisfies 


t 

TA  = exp[i  J dt'  A (t1)  ] 
<1 


(4.11) 


X>$A  = 0 (4.12) 

and  the  normalization  condition  N = (Ya|  Ta)  gives 

t i 

exp[-2  J dt'  {/m(  A (<'))}  ] = ( ($A^))  (4-13) 

h 

Since  Re(A)  does  not  appear  in  the  normalization  constraint  in  Eq.(4.13),  we  are 
free  to  choose  it  equal  to  zero.  Then  the  norm  conserving  solution  to  Eq.(4.5)  is 
given  by 


T a = JV* 


(4.14) 


(4aI*a)s 

where  $a(0  satisfy  Eq.(4.12),  and  the  normalization  constant  N is  constructed 
from  the  initial  conditions. 


4.3.  Time-dependent  Self-consistent  Field  (TDSCF)  Approximation 

Expressing  the  norm  conserving  time-dependent  variational  principle  (TDVP), 
of  the  previous  section,  in  a compact  form  as 


(<TaI»Ta  + ATa)  = 0 (4.9) 

The  frequency  operator  Q,,  which  governs  the  time  evolution  of  the  exact  amplitude 
|Y(i)),  is  of  the  form 
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0 = £ exp 


i 


t 


/ 

V 


= r]  exp 


(4.20a) 

(4.206) 


we  remove  the  normalization  constraint  A from  the  differential  equations  corre- 
sponding to  £'(0  and  p'(t),  respectively.  We  find  the  differential  equations  corre- 
sponding to  the  self-consistent  primary  and  secondary  region  amplitudes  £(t)  and 
rj(i),  respectively,  in  their  final  form  as 


it  - 6pc  - 

ip  — Qsp  — 


(v\*i) 

(C|ftp3IC) 

(CIO 


C = o 


p = o 


(4.21a) 

(4.216) 


Imposing  the  normalization  constraint  N = (T \(t)  | T and  using  the 
transformation  equations  (4.19a,b)  and  (4.20a,b),  we  write 


l 

N = (OHO?)  exp[~  J { 


(v  10  + (v\v) 
(v\n) 


(4.22) 


+ 


(CIO  + (cio 
(CIO 


i 

}dt'  — 4 J dt' Im  {\  (i1))] 


Subtracting  complex  conjugate  of  Eq.(4.18a)  from  Eq.(4.18a),  and  using  the  result 
in  Eq.(4.22),  we  express  the  normalization  constraint  as 
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exp 


i 

J dt'lm  {A  (t1)  } 


N 


(Zvltn) 


(4.23) 


where  again  as  in  the  previous  case  we  are  free  to  choose  Re(A)  = 0,  since  the 
normalization  constraint  does  not  involve  Re(A).  Finally,  using  the  transformations 
in  equations  ( 4.19a, b)  and  (4.20a, b)  in  the  TDSCF  approximation  in  Eq.(4.17), 
we  get 


N = £77  ex 


l 

p[+i  J { 


i(v\v)  - {n\fr\v) 
ii\n) 


(4.24) 


mi)  + gi^io  _ 2Jm  (A  y))}  dt' ) 

m 


where  again  we  use  Eq.(4.18a)  to  write 


Hv\n)  - 
tflv) 

m 


(4.24) 


(tv\npam 


- A 


which  together  with  the  normalization  constraint  in  Eq.(4.23),  and  the  choice  Re(A) 
= 0,  gives  the  norm  conserving  TDSCF  approximation  in  its  final  form  as 


Tj  = jvi  — C.  -!!- -exp 

«l«-  (nW 


f , {(niW’Kn) 

- 1 1 


(4.26) 


The  amplitudes  £(t)  and  r/(t),  corresponding  to  the  self-consistent  primary  and 
secondary  region  dynamics,  respectively,  satisfy  the  following  coupled  equations: 
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ii-freff(t)t  — o 

in  - &eff  (O7?  = 0 


(4.27a) 


(4.276) 


where  the  time  dependence  in  the  effective  frequency  operators  Vtve^(t)  and  (lse^(t) 
are  given  by 


(4.28a) 


(4.286) 


Similar  to  the  time-dependent  Hartree  (TDH)  factorization  of  the  total  wave  function 
[Schatz,  77;  Gerber,  Buch  and  Ratner,  82;  Kugar,  Meyer  and  Cederbaum,  87; 
Makri  and  Miller,  87]  for  a coupled  many-mode  problem,  we  also  achieve  a formal 
separation  of  the  total  amplitude,  T (t),  into  self-consistent  primary  and  secondary 
region  factors  £(t)  and  t]( t),  respectively.  The  time-dependence  in  each  of  these 
factors  is  governed  by  an  effective  frequency  operator,  which  includes  an  average 
of  the  interaction  frequency  operator,  over  the  other  factor.  The  essential  feature  of 
the  TDSCF  approximation  is  that,  while  it  provides  the  formal  factorization  of  the 
dynamics  of  a multivariable  system  into  primary  and  secondary  region  dynamics,  it 
still  permits  the  self-consistent  energy  exchange  between  these  two  regions  through 
the  effective  couplings  in  the  frequency  operators.  An  essential  requirement  in  our 
work  is,  that  we  must  have  a norm  conserving  form  for  the  amplitude  T(r)  for  a 
complex  time  r > 0,  and  a correct  phase  factor  in  Eq.(4.26).  This  is  because  we  are 
interested  in  computing  the  overlap  of  the  complex  time  amplitude  with  its  initial 
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value  at  r = 0;  therefore  the  phase  factor  an  integral  of  the  effective  coupling  from 
0 to  r,  does  not  get  canceled  by  its  own  value  at  r = 0.  The  norm  conserving  form 
for  the  amplitude  in  Eq.(4.26)  is  essential  because  the  effective  frequency  operators 
Clpeff(T)  and  are  non-hermitian  for  any  complex  time  r > 0. 

4.4.  Primary  and  Secondary  Region  Factorization  of  TCFs 


Using  the  TDSCF  factorization  of  the  amplitude  T(r)  as  shown  in  Eqs.(4.26- 
28),  in  Eq.(3.34)  of  the  previous  chapter,  we  write  the  complex  time  molecular 
TCF  as 


^>W  = £»W/;Ms/(t)  (4.29a) 

/ 

where  we  have  replaced  the  operator  A by  a molecular  dipole  operator  D,  restored 
the  subscript  I corresponding  to  the  sum  over  the  initial  nuclear  vibrational  state 
| xpj),  and  the  complex  time  r is  given  by  t = t — ihu.  We  define  the  “primary 
region  TCF”  //(r),  as  a dimensionless  overlap 


hij) 


(t'(0)l(f00) 


(4.296) 


and  the  “secondary  region  TCF”  gi(r),  also  as  a dimensionless  overlap 


_ j7!1  (0)  W (r)) 

WW) 


(4.29c) 


The  dimensionless  primary  region  TCF  //(r)  and  the  secondary  region  TCF  gi(r), 
together  with  an  overall  phase  factor  a/(r),  which  is  an  integral  of  the  self- 
consistent  coupling  from  0 to  r given  by 
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Oil 


(eVl^'W)  , 

— — dt 

<€VICV> 


(4.29) 


o 

completely  describe  the  time-dependence  in  the  total  TCF  F ^ t ).  The  normaliza- 

tion constant  Nj  in  Eq.(4.29a),  which  is  computed  from  the  given  initial  conditions 
£/(0)  and  77/(0)  as  according  to 


Nj  = (e7  (0)  rj1  (0)  |£7  (0)  r/1  (0))  , (4.29d) 

provides  correct  dimensions  to  the  overall  TCF. 


CHAPTER  5 

TIME-DEPENDENT  EFFECTIVE  FREQUENCY  OPERATORS 


In  the  present  chapter,  we  consider  a general  nuclear  vibrational  Hamiltonian 
in  a given  electronic  state,  and  solve  the  equations  of  motion  for  the  secondary  re- 
gion dynamics.  In  section  5.1,  we  describe  the  time-dependent  effective  frequency 
operators  by  considering  many  degrees  of  freedoms  in  the  secondary  region  inter- 
acting strongly  with  a single  degree  of  freedom  in  the  primary  region.  We  model 
the  strongly  coupled  primary  and  secondary  regions  couplings  which  are  of  arbri- 
trary  analytic  form  in  the  variables  of  the  primary  region,  and  have  a linear  and 
bilinear  dependence  in  the  variables  of  the  secondary  region.  Section  5.2  describes 
the  effective  frequency  operators  for  a weak  coupling  limit,  and  ignores  the  terms 
in  the  coupling  which  are  bilinear  in  the  variables  of  the  secondary  region.  The 
solutions  to  the  equations  of  motion  for  the  secondary  region  operators,  in  both 
of  these  cases,  are  reduced  to  analytic  forms,  which  are  then  used  to  construct 
the  propagators  for  the  secondary  region  dynamics.  The  primary  region  dynamics 
will  be  described  in  the  next  chapter  by  using  a numerical  path-integral  method. 
In  section  5.3  we  derive  an  operator  equation  of  motion  for  the  primary  region 
dynamics,  show  its  relationship  to  the  well  known  phenomenological  Generalized 
Langevin  equation  (GLE)  [Mori,  65],  by  constructing  the  fluctuation  force  term  and 
the  dissipative  memory  kernel,  and  prove  the  Fluctuation-Dissipation  Relationship 
(FDR)  [Van  Hove,  54]  explicitly  for  a general  intermolecular  process.  Section  5.4 
considers  many  degrees  of  freedom  in  the  primary  region  dynamics  interacting  with 
a set  of  coupled  harmonic  oscillators  in  the  secondary  region  dynamics. 
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5.1.  Intramolecular  Processes:  Strongly  Coupled  Primary  and  Secondary  Regions 

In  a general  intramolecular  photointeraction  process,  following  the  absorption 
of  light  in  the  visible  and  UV  region,  we  consider  electronic  transitions  between 
fixed  electronic  states  i and  j.  Therefore,  for  notational  brevity  throughout  this 
chapter,  we  avoid  writing  i and  j as  the  specific  references  to  the  initial  and  the  final 
electronic  states,  respectively.  Consider  a general  nuclear  vibrational  Hamiltonian 
H in  the  given  electronically  excited  state  as  consisting  of  the  following  terms: 


p2  , , 

Hv  = -^-  + X°  (X) 
2m  x V / 


(5.1a) 


is  the  Hamiltonian  for  the  primary  region  dynamics  with  an  arbitrary  potential 
I/0(X)  given  either  in  a parametrized  analytic  form  based  on  experimental  infor- 
mation, or  on  a grid  in  coordinate  space  from  ab  initio  or  semi-empirical  electronic 
structure  calculations; 


Pi 


1 


= £ itk  + 


(5-16) 


is  the  Hamiltonian  for  the  secondary  region  expressed  as  a sum  of  k(k  = 1, 3...N) 
harmonic  oscillators,  such  that  the  mass  m*  and  the  frequency  w*  of  the  k1h 
harmonic  oscillator  correspond  to  the  mass  and  the  frequency  of  the  kth  normal 
mode  vibration  in  the  molecule.  The  coupling  between  the  primary  and  the 
secondary  regions  is  expressed  as 


H”’  = Y,  { PkhxPx  + Vk  (x)  <?*  + £ Qt'Vl'k  (x)  Qt ] 
k l k'  ) 


(5.1c) 
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where  the  arbitrary  functions  Vk(X)  and  Vkk'(X)  are  operator  dependent  forces 
and  potential  curvatures,  respectively,  and  the  bilinear  momentum  coupling  pa- 
rameter Xkx  has  dimensions  of  inverse  mass.  The  arbitrary  functions  Vk(X ) and 
Vkk'(X ) are  also  given  either  in  an  experimentally  parametrized  analytic  form,  or 
on  a grid  in  coordinate  space  as  computed  in  an  electronic  structure  calculation. 
Following  the  absorption  of  light  in  the  visible  or  UV  region,  one  finds  that  the  pri- 
mary region  dynamics,  describing  for  example  rearrangement  in  the  system  due  to 
the  formation  or  breaking  of  a bond,  is  strongly  coupled  to  the  rest  of  the  molecule. 
Therefore,  we  have  added  to  the  coupling  the  terms  linear  and  bilinear  in  the  oper- 
ators of  the  secondary  region.  We  also  note  that  the  formulation  does  not  depend 
upon  the  strength  of  the  various  terms  in  the  coupling,  and  in  the  next  section  we 
consider  the  weak  coupling  limit  by  ignoring  the  first  and  the  last  terms  in  Eq.(5.1c). 

For  a single  degree  of  freedom  in  the  primary  region  dynamics  coupled  to  many 
degrees  of  freedom  in  the  secondary  region,  we  write  the  norm  conserving  TDSCF 
approximation  for  the  amplitude  T 7(tf)  as 


T7  (t)  = Nf  e+,a/W 


£7(0 


N 

n 


•»*'(*) 


(£7I£7)’  til  ( rikl\rihI)2 


(5.2) 


Here  the  secondary  region  amplitude  77 7(f)  in  Eq.(4.26)  of  the  previous  chapter  is 
written  as  a product  over  k,  and  the  amplitude  rjkl(t ) corresponds  to  the  kth  normal 
mode  motion  in  the  secondary  region  dynamics.  The  time-dependent  effective 
frequency  operator  for  the  primary  region  dynamics  is  given  by 
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fr'ff  (0  = l{Hp  - &,  + E<P*  (i,  hi)  AtA-PA- 

* 

+ V‘(A)  qt  (i,  [17])  + iv“(X)  si  (t,  hi)  (5.3) 

5 E ?*'  <<>  M)  V*  *(A-)g*  (t,  hi)  - <}] 

^ k^k1 

where  we  have  omitted  specific  indices  within  the  functional  dependence  for 
notational  brevity,  and  on  the  ith  potential  energy  surface  we  have  written  the 
stationary  state  eigenvalue  Ej  as  Ej  = Ep  + Ej  + Eps , with  Ej=  ^2kEj, 
and  Eps  = Yhk  The  time-dependent  coefficients,  which  couple  the  primary 
region  dynamics  to  the  normal  mode  motions  in  the  secondary  region,  are  the 
normalized  expectation  values  of  the  corresponding  operators  within  the  amplitudes 
coresponding  to  the  secondary  region  dynamics.  For  example,  the  coefficient 
ql(t,  [?/])  is  of  the  form 


«i(*.fo]) 


h^lQlh*7) 


(5.4) 


The  time-dependent  effective  frequency  operator  Clse^(t)  for  the  motion  in  the 
secondary  region  dynamics  is  given  by 

6!//«  = E“«//W  <5-5) 

k 

where  the  effective  frequency  operator  for  the  kth  normal  mode  motion  in  the 
secondary  region  is  written  as 
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&cf/(t)  = + ^mivl(t,[(])Ql  + Pt\kxpx 

+ ®*  (<,  K,  «))  <?*  + £*  (*,  K,  5])  - £7  - £?*]  (5.6a) 

where  77  means  any  77*/  with  fc'  ^ fc.  The  time-dependent  effective  frequency 
[£]),  the  effective  force  vk(t,  [£,  77] ),  and  the  effective  energy  £k(t,  [£,77]) 
between  the  kth  normal  mode  in  the  harmonic  bath,  the  primary  region,  and  the 
rest  of  secondary  region,  are  written  as 

^(Me])={^+”tt^’K1)};  (5.66) 

®‘  («,  KD  = ■>*  (Mf])  + \ E I"*'*  ('•  1®  + (*-  ({])}  «•  («,  M) 

z k'jtk 

(5.6  c) 


£k  (*,  [£, »/])  = X!  (*»  M)  (<,  [£]) 

k'^k 

+ ^ (*.  [£])  ?Jb'  (*,  [»?])  + |*>*'*'  (f,  [£])  q2k,  (f,  [77])  (5.6c?) 

+ \ qk"  (<»[^])t,*"fc,(<»K])9fc'(f>  W) 

Z k"^k'jik 


57 


The  time-dependent  coefficients  that  couple  the  kth  mode  of  the  secondary 
region  to  the  primary  region  dynamics  are  defined,  for  example,  by 


»*'*  (<,(«])  = 


(x)  |€J> 


(5.6e) 


The  time-dependent  effective  phase  a(t,  [r/,^]),  defined  in  Eq.(  4.29d)  and  needed 
to  compute  the  total  TCF  F ^ ^(t)  in  Eq.(4.29a),  is  given  by 

x 

(t,  [£,  */])  = jr  J dt'  (<',  [*?])  hxPX  (*',  [£]) 

0 k 

+ »*  (<’,  Kl)  it  (*',  W)  + (<’,  K])  ?1  («,  W)  (5.7) 

+ 5X>  (AfilK  V.  lfl) « (<.  M)} 

The  form  of  the  effective  frequency  operator  hkejj(t)  in  Eq.(5.6)  for  the  kth  mode 
of  the  secondary  region  dynamics,  is  that  of  a linearly  forced  parametric  oscillator. 
Operators  Qk(t)  and  Pk(t),  corresponding  to  the  kth  mode  motion  in  the  secondary 
region  dynamics,  evolve  according  to  the  following  equations  of  motion 


Qk  = * 
Pk  = i 


a 

(t) , Pk 


efftt)  iQk 

k 

eff 


(5.8a) 
(5.8  b) 


where  [...]c  are  quantum  commutators.  Using  Clk^(t)  from  Eq.(5.6a),  it  is  easy 
to  show  that  Qk(t)  satisfies 
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Qk  + uk  (t)  Qk 


*kXPX  (t)  - vk  (t) 
mk 


(5.9) 


where  we  have  dropped  the  functional  dependence  on  £ and  fj  for  notational  con- 
venience. Equation  (5.9)  is  the  equation  of  motion  of  a linearly  forced  parametric 
oscillator,  therefore  we  solve  it  by  separating  it  into  a homogeneous  equation  and 
a nonlinear  forcing  term.  We  are  using  this  approach,  because  it  is  suitable  for 
application  to  computations  in  the  complex  time  plane,  and  also  because  it  is  use- 
ful in  deriving  a parallel  between  the  TDSCF  dynamics  and  the  GLE  approach  to 
stochastic  dynamics. 

We  introduce  uk(t)  and  i ’k(t),  two  linearly  independent  solutions  of  the  ho- 
mogeneous equation 


d2  2 

dti+Uk(t) 


Uk(t 
vk  (t. 


= 0 


(5.10) 


with  the  initial  conditions  uk(0)  = 0,  and  u*(0)  = 1.  The  full  solution  to  the 
operator  Eq.(5.9)  is  written  as 


Qk  (t)  = Qk  (0)  vh  (<)  + Pk  — uh  (<)  + yk  ( t ) (5.11a) 

mku>k 

where  Qit(0)  and  Pk( 0)  are,  respectively,  the  initial  values  of  the  operators  Qk(t) 
and  Pk(i)  at  the  initial  time  t = 0,  and  the  the  last  term  yk(t ) is  an  integral  of 
the  form 


yk  (t)  = ~^~k  J dt'  {^itXPAT  (0  - vk  (*')}  uk  (t  - t')  (5.116) 

o 
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Using  equations  (5.11a)  and  (5.11b)  with  the  definitions  in  Eq.(5.4)  for  the  time- 
dependent  coffecients  of  the  coupling  between  the  primary  and  the  secondary  region 
dynamics,  we  express  the  coefficients  needed  to  compute  the  effective  frequency 
operator  f)^-y(t)  for  the  primary  region  dynamics  as 

Qk  (0  = Qk  (0)  vk  (t)  + EkifU ujt  (t)  + yk  (t)  (5.12a) 

mkiOk 


Pk  (t)  = Pk  (0)  uk  (t)  + mhqk  (0)  Vk  (*)  + rnkVk  (t)  ~ hx™kPX  (<) 

(5.126) 


Qk  (*)  — qI  (o)  vl  (0  + Pk2  \ul  (0  + yl  (0 

mkuk 


+ 2 [qk  (0)  vk  ( t ) + (<)]•  Vk  (t)  (5.12c) 

m.LUJL 


where  the  initial  values  y*(0),  p*(0),  and  g£(0)  are  defined  as  according  to 


a,-,  <■?*'(<>)  Iff  totJ«>)) 
qti>  (Vl,(0)\rik,m 


(5.12  d) 


At  any  time  t the  time-dependent  effective  frequency  operator  £lpefj(t)  depends 
upon  the  previous  history  of  the  primary  and  the  secondary  region  dynamics 
through  the  term  y*(<),  which  involves  integral  of  the  nonlinear  forcing  term  from 
0 to  t,  therefore,  in  our  TDSCF  approach  we  refer  to  yjt(f)  as  a “memory”  term. 
Furthermore,  in  a later  section  of  this  chapter,  we  show  a direct  relationship  between 
the  TDSCF  approach  and  the  GLE  approach  to  stochastic  dynamics,  by  deriving  a 
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fluctuation-dissipation  relationship  between  the  memory  kernel  and  the  fluctuation 
force. 

It  appears  from  Eq.(5.11b),  that  the  computation  of  the  memory  term  yk(t ) 
at  any  time  t,  involves  knowing  the  integrand  at  all  previous  times  t'  < t,  as 
well  as  at  the  present  time  t’  = t.  It  is,  however,  not  difficult  to  see  that  the 
choice  of  the  initial  conditions,  uk(0)  = 0 and  u*(0)  = 1,  for  the  solutions  of  the 
homogeneous  Eq.(5.10),  make  the  integrand  in  the  memory  term  yk(t)  equal  to  zero 
at  t'  = t.  Therefore,  this  choice  of  the  initial  conditions  is  significant  in  the  sense 
that  it  enables  us  to  perform  the  TDSCF  computation,  for  a nonlinearly  coupled 
primary  and  the  secondary  regions  dynamics,  without  going  through  iterations  at 
each  time.  In  other  words,  these  initial  conditions  have  made  it  possible  to  perform 
our  computations  in  a step  by  step  progression  along  a desired  axis  in  the  complex 
time  plane,  while  at  the  same  time  maintaining  the  self-consistency  in  the  dynamics. 

We  also  note  that  the  linearly  independent  solutions  uk(t)  and  vk(t)  which 
satisfy  the  homogenous  Eq.(5.10)  for  the  kth  normal  mode  motions  in  the  secondary 
region,  are  completely  independent  of  the  simultaneous  solutions  uk>(t)  and  vk>(t ) 
for  all  k'  / k.  Therefore,  for  any  secondary  region  made  up  of  N coupled  harmonic 
oscillators,  we  have  to  simultaneously  solve  only  2 N uncoupled  homogeneous 
equations  of  motion  to  achieve  the  self-consistent  solution  of  the  full  dynamics. 

The  initial  values  qk( 0),  and  pk(0)  defined  as  according  to  Eq.(5.12d),  are  the 
normalized  expectation  values  of  the  corresponding  operators  within  the  amplitude 
rjkI(0).  They  can  also  be  thought  of  as  the  mean  position  and  momentum  cor- 
responding to  the  kth  normal  mode  motion  in  the  secondary  region  dynamics.  If 
the  secondary  region  is  initially  prepared  at  thermal  equilibrium  in  the  presence  of 
the  primary  region,  then  we  can  use  a suitable  thermal  distribution  for  these  ini- 
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rial  values  to  introduce  a temperature  dependence  in  the  dynamics  of  the  primary 
region  operators. 

We  end  this  section  by  using  equations  (5. 5-5. 6)  for  the  definition  of  the  effec- 
tive frequency  operator  Cl*jj(t)  for  the  kth  mode  motions  in  the  secondary  region 
dynamics,  and  equations  (5.11-5.12)  as  the  solutions  for  the  operator  Eq.(5.9),  to 
construct  the  propagator  for  the  kth  mode  motion  in  the  secondary  region  dynam- 
ics [Yamashita  and  Miller,  85]  as 


(<?*,&, 


I r 

ihuk(t)J  P 2 huk(t) 

t 

{(Ql  + QtH  (<)  - 2 QkQ'k  - f h ((')  ut  ((')  dt' 

mkuk  J 
0 
t 

— J fk  ( t ')  uk  (t  - t')  dt'  (5.13a) 


mkUk 


mWk 


l L 

J fk  ( t ')  uk  ( t - t')  dt'  J fk  (t")  uk  (t")  dt"  } 
0 0 

+ L(E$  + E*k)t] 


where  the  effective  force  /*(<),  in  the  integrand,  is  given  by 


fk(t)  = mk\kXpx(t)-vk(t ) (5.136) 

and  the  time  evolution  of  the  amplitude  corresponding  to  the  kih  normal  mode 
motion  in  the  secondary  region  dynamics  is  described  by 

+oo 

J dQ'k  (Qk\Kk  (Qk,Q'k,t)  \Q'k)  r,k>  (Q'k,0) 

— oo 


(5.14) 
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This  is  about  as  far  as  one  can  proceed  in  general,  i.e,  without  numerically  solving 
the  homogeneous  Eq.(5.10)  for  a given  process. 


5.2.  Intermolecular  Processes:  Weak  Coupling  Limit 


Intermolecular  photointeraction  processes  such  as  overtone  line  broadening, 
dephasing,  and  energy  transfer  in  an  extended  molecular  system,  following  a 
localized  absorption  of  light  in  the  visible  and  UV  region,  are  often  described 
by  the  time  evolution  of  a localized  primary  region,  interacting  weakly  with  a 
collection  of  coupled  harmonic  normal  mode  motions  in  the  secondary  region.  In 
the  formulation  of  the  previous  section,  we  drop  the  bilinear  coupling  terms,  and 
keep  the  term  which  is  linear  in  the  operator  dependence  of  the  secondary  region 
dynamics  and  has  an  arbitrary  analytic  form  in  the  operator  dependence  of  the 
primary  region.  In  other  words,  the  effective  frequency  operator  &peff(t)  for  the 
primary  region  dynamics  simplifies  to 


Hp- 


£?+£{v‘(*) 

k 


(5.15) 


where  the  effective  coupling  of  the  primary  to  the  secondary  region  is  also  defined 
by  the  example  in  Eq.(5.4).  The  time-dependent  effective  frequency  operator  for 
the  kth  mode  of  the  secondary  region  simplifies  to 


21  + I mwlQl  + vk  «,  [£))  Qk  - (e?  - Ef) 


(5.16) 
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and  the  time-dependent  effective  phase  needed  to  compute  the  total  TCF  is  written 
as 


t 

«(<»[£>*/])  = J dt>  {«*(<>[£])  9*  (M»d)  ~Eik}  (5-17) 

k o 

The  form  of  the  operator  in  Eq.(5.16)  is  that  of  a linearly  forced  harmonic  oscillator. 
It  is  not  very  difficult  to  see  that  the  homogeneous  Eq.(5.10)  for  the  functions 
Uk(t),  and  t’jt(f)  reduces  to  the  equation  of  motion  for  a harmonic  oscillator  with 
frequency  u;*  and  mass  m*.  Therefore,  the  functions  uk(t)  and  u*(f)  with  the 
initial  conditions,  u*(0)  = 0 and  u*(0)  = 1,  simply  reduce  to  u*(£)  = sin(u;fcf)  and 
Vk{t)=  cos (0;^#),  respectively.  The  time-dependent  effective  coupling  <?*(*,  [77]), 
which  couples  the  primary  region  dynamics  to  the  kth  mode  motions  in  the 
secondary  region  dynamics  is  given  by 

qk(t)  = 9Jb(0)coswfc<  + — — ■ sin ukt  + yk  (t)  (5.18) 

rrikUk 

where  the  “memory”  term  yk(t)  for  the  weak  coupling  case  simplifies  to 

t 

yk  (t)  = f dtt'  vk  (t1)  sinwjt  (t  — t')  (5.19) 

rwjfcWfc  J 

0 

Here  again,  we  notice  that  the  computations  for  the  memory  term  in  Eq.(5.19) 
involve  the  time  dependence  in  the  effective  potential  vk(t)  only  for  the  previous 
times  t'  < t.  Therefore,  the  TDSCF  computations  for  the  primary  region  dynamics 
can  be  performed  in  a step  by  step  progression  along  any  chosen  time  axis. 

Finally,  the  propagator  for  the  secondary  region  dynamics,  in  the  weak  coupling 
case,  is  expressed  in  an  analytic  form  as  [Schulman,  81] 
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Kk  (Qk,Q'k,i)  = 


( ) 3 „,,[<?!*_{  (Ql  + ft)  cos utt 

\2nih  sinwjtf / y l2h  sinu,’*f  1 \k  k J 

/ 2Qa-  f jj  k u\  . ±<  ^Qk 


- 2Q*g'*  - 


sin  — 


mk^k 


mWk 


(5.20) 


J dt'vk  (t')  si 

o 

t 

x J dt'vk  (t1)  sincjfc  (<  — — 

o 

t t' 

x J dt'vk  (t1)  sinu;*;  ( t — t')  J dt"vk  ( t ")  sin  } 

o o 

+ l (£/  + Ef ) t } 


5.2.  Relationship  to  the  Generalized  Langevin  Equation  Approach 


The  method  of  factoring  the  dynamics  of  a large  molecular  system  into  a 
“primary  region”  consisting  of  a few  important  degrees  of  freedom,  interacting 
with  the  rest  of  the  system  labeled  as  the  “secondary  region”,  is  well  known  in 
statistical  mechanics,  where  one  often  needs  to  deal  with  the  dynamics  of  a “system” 
in  thermal  equilibrium  with  a large  heat  reservoir  called  “bath”.  Among  many 
different  approaches  to  deal  with  the  factorization  of  a large  system  into  coupled 
primary  and  secondary  region  dynamics,  the  approach  based  on  the  Generalised 
Langevin  Equation  (GLE)  [Mori,  1965],  solving  it  under  different  approximations 
for  the  primary  region  dynamics,  is  well  covered  in  the  literature.  The  terms 
which  affect  the  dynamics  of  primary  region  variable,  due  to  the  coupling  with  a 
large  heat  reservoir,  are  the  “fluctuation  force”  term  and  the  “memory”  term.  The 
form  of  the  fluctuation  force  term  is  decided  by  the  dynamical  and  the  statistical 
properties  of  the  secondary  region,  and  depends  largely  upon  the  way  the  secondary 
region  is  brought  to  a thermal  equilibrium  in  the  presence  of  the  primary  region. 
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If  the  primary  region  is  coupled  to  an  infinite  heat  reservoir  for  a sufficiently  long 
time,  the  primary  region  also  attains  a steady  state  due  to  a relationship  between 
the  memory  kernel  and  the  flutuation  force,  this  relationship  is  often  referred  to 
as  the  Fluctuation-Dissipation  Theorem  [Van  Hove,  54;  Berne,  71].  The  GLE 
is  a stochastic  differential  equation  with  additive  fluctuation  force  and  dissipative 
memory  terms. 

Starting  with  the  time-dependent  effective  frequency  operator  £tpeff{t)  in 
Eq.(5.15),  we  derive  an  operator  equation  for  the  primary  region  dynamics,  show 
its  relationship  with  the  GLE  by  constructing  the  fluctuation  force  and  the  memory 
kernel  terms,  and  finally  derive  the  Fluctuation-Dissipation  Relationship(FDR)  for 
our  TDSCF  dynamics. 

The  time  evolution  of  the  operators  X and  Px  for  the  primary  region  dynamics 
is  described  by  the  following  operator  equations 


X = i 

Px  = i 


npeff(t),x 
npeff(t),P x 


(5.21  a) 
(5.21  b) 


Using  Eq.(5.15)  for  the  effective  frequency  operator,  Eq.(5.21b)  is  written  as 


= i [y0  (*) -rir]c+£i  [v*  (*)  • PA « «* « <5-22> 


where  we  avoid  the  functional  dependence  in  q^t)  for  notational  simplicity.  From 
Eq.(5.18)  the  effective  coupling  qk(t ) is  written  as 
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Qk  (t)  = I qk  (0)  + — - — vk  (0)  ) cos  wkt  + 


mhuj  k 


Pk  (0)  . 

1 sin  ujkt 

rrikUk 


vk(t ) 1 


mkUk  m*u;£ 


l 

J dt'vk  (<')  cosa )k  (t  - t')  (5.23) 


where  we  have  used  integration  by  parts  in  the  memory  term  in  Eq.(5.19),  and 
substituted  the  result  in  Eq.(5.18).  Using  Eq.(5.23)  in  Eq.(5.22)  and  rearranging 
the  result  we  get 


Px -<  [n""J (x,Px,t),Px]c  = ^ [v* (*) ’Px 

rJ{vkP)^i) 


X J dt' 
0 


2rrikmxu>l 

+ Rk  (<)  } 


COS  UJk  (t  — t') 

(5.24) 


where  the  commutator  involving  the  modified  frequency  operator  Cl^od(t)  given  by 


»u°d  (*M  - 1 + v*  (*)  - E Vi3vt  (*)  - E ■ 


PI 


vk(t) 


2mx 


mku  % 


(5.25  a) 


describes  the  time-dependence  in  the  uncoupled  primary  region  dynamics  in  the 
presence  of  a stationary  secondary  region,  and  the  fluctuation  force  Rk(t)  is  given 
by 


Rk  (<)  = jr  L 


V*  (x)  ,Pv]c  | (« (0)  + COSU ,kt  + sin utt  } . 


mkuJk 


mkujk 


(5.25  b) 
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The  integrand  in  Eq.(5.23)  has  been  simplified  to  the  form  in  Eq.(5.24)  by  using 
the  equation  of  motion  of  the  operator  Vk(X)  to  write 


The  form  of  Eq.(5.24),  a general  operator  equation  for  the  primary  region  dynamics 
in  the  presence  of  a weakly  coupled  secondary  region,  is  an  analogue  of  the  GLE 
in  our  TDSCF  approach.  We  have  considered  a coupling  which  is  linear  in  the 
variables  of  the  secondary  region,  and  has  a general  analytic  form  V*(X)  in 
the  variables  of  the  primary  region.  Therefore,  Eq.(5.24)  for  the  primary  region 
dynamics  in  our  TDSCF  approach  is  more  general  than  the  usual  GLE  given  in  the 
literature.  As  an  example,  we  consider  a simple  case  of  Vk(X ) = r*Ar,  with  Tk 
a free  force  constant,  and  use  this  in  Eq.(5.24)  to  write  it  as 


(5.25c) 


with 


/f  i \ _ (e(Q)i[-]cie(Q)) 
u"‘u  <£(o)iao)> 


(5.25  d) 


Px  - i [£lMod  (X,  Px,  t)  , Px\  c - J dt'M  (*-*')(  ) = *(*) 


o 


(5.26) 


where 


(5.27a) 


k 


and 
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R(t-t')  = '£,Rk(t~t') 


(5.27  b) 


k 


with 


mkuk 


(5.28a) 


and 


x (0)  ) cos  ( uiki ) + 


respectively.  Equation  (5.26)  is  the  analogue  of  the  usual  GLE  for  a bilinearly  cou- 
pled system-bath  problem,  where  the  fluctuation  force  term  R(t)  and  the  memory 
kernel  — are  given  in  Eqs(5.27a)  and  (5.27b),  respectively.  It  is  not  difficult 
to  see  that  the  primary  region  dynamics  depends  upon  the  initial  values  of  the  sec- 


R{t), which  contains  the  normalized  expectation  values  ?*(())  and  pk( 0),  given  for 
example  by 


Therefore,  the  fluctuation  force  term  R(t)  can  be  given  a stochastic  interpretation 
by  observing  that  these  expectation  values  are  described  by  a distribution  function 
which  depends  upon  the  initial  temperature  of  the  system.  If  the  initial  state  of 
the  secondary  region  is  an  equilibrium  state,  then  the  distribution  of  all  the  bath 
operators  is  the  same  and  consequently  R(t)  is  Gaussian  (Central  Limit  Theorem). 
The  mean  value  of  the  resulting  Gaussian  depends  upon  the  form  of  the  distribution 
function.  If  initially  the  secondary  region  has  been  brought  to  the  equilibrium  in  the 


ondary  region  operators  Qjt(0)  and  Pk( 0)  only  through  the  fluctuation  force  term 


(5.29) 
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presence  of  the  primary  region,  a suitable  choice  of  the  distribution  function,  which 
give  the  mean  value  of  the  fluctuation  force  equal  to  zero,  is  a displaced  Boltzman 
distribution  function  [Zwanzig,  73;  Lindenberg  and  Seshadri,  82]  expressed  as 

IV(q(0),p(o)),/3)  = ^exp(-/3^[  ^ + Vw|{?t(o)  + ^x(o)}’]) 

(5.30) 

The  thermal  averaged  correlation  function  of  the  force  term  R(t)  then  is  written  as 

(R(t)R(t%„  = E E / ‘W  f dp(°)  w(q(0),  P(0)),  /?) 

J J (5.31) 

x Rk(qk(0),Pk(0),t)Rk'(qk'(0),Pk'(0),t) 

where  -Rfc(5*(0)>Pit(0),f)  and  Rk'(qk'(Q),Pk'(Q),t')  are  as  defined  in  Eq.(5.28b). 
Using  the  distribution  function  W(q(0),  p(0))  in  Eq.(5.31),  we  obtain  the  correla- 
tion function  of  the  fluctuation  force 


<* ««((')),„  = kBT E c°»u>*  (t  - 1‘) 

k mk^k 


which  is  the  same  as 


(5.32) 


(. R(t)R(t'))'U  = kBTM(t-t ')  . (5.33) 

This  shows  a FDR  between  the  memory  kernel  M(t  — t1)  and  the  fluctuation 
force  R(t).  For  the  more  general  operator  equation  in  (5.24)  for  the  primary 
region  dynamics,  we  get  a similar  relationship  between  the  memory  kernel  and  the 
fluctuations  in  the  force. 
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5.4.  Treatment  of  Many  Degrees  of  Freedom  in  the  Primary  Region  Dynamics 


To  treat  many  degrees  of  freedom  l = 1 toL  in  the  primary  region  dynamics, 
we  generalize  the  Hamiltonian  operators  of  section  5.1,  and  write 

= E ||  + V°&  ■ (5.34a) 


/=! 


A‘  = E it + lmk“kQl 


h=l 


2m  ^ 2 


(5.336) 


H”  = E<E  P^‘P'  + VtW&  - ? E 0k'vk'k(S)0k  } (5.33c) 

h l L k'^k 

where  X=  {X\,X2 are  the  operators  of  position  variables  of  the  primary 
region,  and  k runs  over  all  the  harmonic  modes  in  the  secondary  region  dynamics. 
Using  the  TDSCF  approximation 


T' = 7S77TT II  T^TT^r  «p(*“(')) 


(5.34) 


(('  I (')i  V (riu  i ’iu'rj 

the  time-dependent  frequency  operator  for  the  primary  region  dynamics  is  given  by 

«?//(*)  = \[H”  - E]  + 


k l 


+ s»(f)v‘(S)  + §4(t)v“(X) 

I E »-(0V*'*(*)qk(t)  - ^ ) ] 


(5.35) 


and  the  time-dependent  effective  frequency  operator  for  the  kth  mode  of  the 
secondary  region  dynamics  is  given  by 


&eff  (0  — 7^^  + ^mk^k  (0  Ql  + ^2  [Pk^klPl  (<)} 

+ vk(t)Qk  + £k(t)-E}-EkjP  ) 


(5.36a) 
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with 


w l + 


mk  . 


(5.366) 


r*  (<)  = (<)  + \Y,  {”*'*  (0  + ykk'  (<)}  <?*'  (0  (5.36c) 


^w  = E p*'  (*)  (o + (o  9t'  (o 

*v*  / 

\ E 9*"  (0”*"*' (*)«'(*)  ) (5.36c?) 

z k"^k'^k 

where  the  time-dependent  coefficients  of  the  coupling  between  the  primary  and  the 
secondary  region  dynamics  are  defined  as  in  sections  5.1  and  5.2.  The  effective 
frequency  operator  for  the  kth  mode  of  the  secondary  region  dynamics  is  still  similar 
to  that  of  an  operator  for  a linearly  forced  parametric  oscillator.  Therefore,  the 
solution  to  the  equations  of  motions  for  the  secondary  region  dynamics,  obtained 
as  described  in  section  5.1  and  substituted  in  Eq.(5.36),  give  the  multi-variable 
time-dependent  frequency  operator  VtpfAt).  One  can  use  Eq.  (5.35)  to  deal  with 
the  dynamics  of  a multi-variable  primary  region  coupled  to  a harmonic  secondary 
region.  If,  however,  the  computations  of  a complex  time  propagator  in  two  or  three 
dimensions  are  prohibitive,  a further  simplification  is  achieved  by  using  an  SCF 
assumption  for  the  primary  region,  so  that 

( rr  e"(.Vi) 

(€*  I €')*  V 


(5.37) 
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in  Eq.(5.34)  to  write  the  TDSCF  approximation  as 

tilt  V\  Jfc/ 


y<  = Nj  TT  J_™_ rr  V(Qk)  (iQ(()) 

' !«")>¥  ■ju>> 


(5.38) 


The  time-dependent  effective  frequency  operator  for  the  primary  region  dynamics 
is  then  written  as 


i 


(5.39a) 


with 


si<//  W - 2^  + «°  (“^'’f)  + St 


+ »*  (*,,«) » (<) + i»“  (*,,«)  ,t2  (<) 
+ ^E«'W'‘,‘(i,’‘)  «« 

z fc'/Jt 


+ (*)  (0  - £/*  } — Ej  ] 


where  v°(Xi,t),vk(Xi,t)  and  vkh(Xi,t ) are  defined  for  example  by 


k.k(  + ^ _ (■..^/-1)/e(,+1)/- 1 vk‘k(k)  i ...e(1-1)Ie(1+1)I-) 


(5.396) 


(5.39c) 


The  effective  frequency  operator  for  the  motion  along  the  1th  degree  of  freedom 
of  the  primary  region  can  be  used  to  compute  its  complex  time  propagator,  and 
the  computational  efforts  for  the  full  primary  region  propagator  scales  linearly  with 
the  number  of  degrees  of  freedom  in  the  primary  region  dynamics,  instead  of 
quadratically,  as  would  happen  if  no  SCF  assumption  was  made  for  the  primary 
region  dynamics. 


CHAPTER  6 

PATH-INTEGRAL  APPROACH  TO  THE  PRIMARY  REGION  DYNAMICS 


The  path-integral  approach  to  the  primary  region  dynamics  is  described  mainly 
by  considering  the  time-evolution  of  the  primary  region  amplitude  |£(f))  according 
to 


i§i\Ht))  = np(t)\((t))  , (6.i) 

where  the  time-dependent  effective  frequency  operator  flp(t)  has  been  described 
in  the  previous  chapter.  We  write  the  integral  solution  of  Eq.(6.1)  in  coordinate 
space  as 


+oo 

(*l£(0>=  J dX0Kp(X,Xo,t)  (Xo|e(0))  (6.2) 

— OO 

where  Kp(X,Xo,t),  the  real  time  propagator  for  the  primary  region  dynamics,  is 
written  as 


+ 00  +oo 

Kp(X,X0,t)=  J cLYjv-i-..  J dXx  (x 


_op 
e N- 


,A  t 


X 


N—l  ^ 


(6.3) 


In  Eq.(6.3),  we  have  divided  the  total  time  t into  N short  time  steps  of  length 
A t =jt  where  N — ► oo.  The  time-dependence  in  the  frequency  operator  Qpt  during 
the  ith  step  in  time  for  i going  from  0 to  N — 1,  is  described  in  Eqs(5.3)  and 
(5.15)  of  the  previous  chapter.  In  the  Lagrangian  formulation  of  the  propagator  the 
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multidimensional  integral  in  Eq.(6.3)  is  also  written  as  a functional  integral,  which 
is  equivalent  to  a summation  over  all  possible  paths  between  the  initial  and  the  final 
positions  in  coordinate  space.  The  number  of  all  such  possible  paths  is  infinite; 
therefore,  since  their  inception  in  the  early  sixties,  the  path-integral  (PI)  methods 
[Feynmann  and  Hibbs,  65],  though  very  appealing  to  the  physical  intuition,  have 
mostly  been  used  as  formal  devices.  In  recent  years,  these  methods  have  also  been 
shown  to  be  a convenient  numerical  tool  in  the  studies  of  many-body  quantum 
mechanical  systems,  and  have  been  used  in  the  computations  of  imaginary  as  well 
as  complex  time  propagators.  All  of  the  recently  developed  numerical  approaches 
to  the  multidimensional  integral  in  Eq.(6.3),  or  equivalently  to  a functional  integral 
in  the  Lagrangian  formalism,  are  in  general  separated  into  two  broad  groups. 
Firstly,  there  are  the  discretized  PI  formulations,  in  the  cartesian  coordinate  space 
[Thirumalai,  Bruskin  and  Beme,  83],  and  in  the  discretized  phase  space  by  using 
the  Fast  Fourier  Transform  (FFT)  algorithms  [Kosloff  and  Kosloff,  83;  Hellsing, 
Nitzan  and  Metiu,  86].  Secondly,  there  are  the  Monte  Carlo  based  approaches  to  the 
multidimensional  integrals;  by  Fourier  Coefficient  path-integration  (FPI)  methods 
[Coalson,  Freeman  and  Doll,  86],  by  random  walk  methods  [Miller,  Schwartz  and 
Tromp;  83],  and  by  using  Metropolis  Monte  Carlo  schemes  [Behrman,  Jongeward 
and  Wolynes,  83]. 

In  section  6.1,  we  briefly  review  a discretized  PI  formulation  in  the  cartesian 
coordinate  space  known  as  the  Numerical  Matrix  Multiplication  (NMM)  method 
[Thirumalai,  Bruskin  and  Beme,  83],  and  discuss  its  shortcomings  due  to  which  it 
becomes  impractical  for  use  in  physical  problems  involving  more  than  one  degree 
of  freedom.  Section  6.2  describes  the  evaluation  of  the  complex  time  memory 
term  y(r)  introduced  in  the  previous  chapter,  and  introduces  a restricted  matrix 
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multiplication  (RMM)  method  to  efficiently  compute  a complex  time  propagator 
for  the  primary  region  dynamics.  By  deforming  the  integration  contour  in  the 
complex  time  plane,  we  describe  in  section  6.3  an  efficient  transformation  from 
real  to  complex  time,  and  show  that  the  efficiency  of  the  RMM  method  can  further 
be  increased  by  using  this  new  transformation. 

6.1.  Numerical  Matrix  Multiplication  (NMM)  Method  and  its  Shortcomings 

For  simplicity,  we  consider  only  a Cartesian  Hamiltonians  in  one  dimension, 
i.e,  the  form  of  the  general  Hamiltonian  is  expressed  as 


are  the  associated  position  and  momentum  operators.  The  solution  to  the  time- 
dependent  Schrodinger  equation  for  the  time-independent  Hamiltonian  of  Eq.(6.4) 
is  obtained  by  writing  the  coordinate  space  propagator  as 


where  [according  to  Thirumalai,  Bruskin  and  Beme;  83,  and  also  for  simplicity]  we 
are  considering  the  propagator  along  the  imaginary  time  axis  such  that  t = — 


(6.4) 


where  mx  is  the  reduced  mass  for  the  degree  of  freedom  along  which  X and  Px 


(6.5) 


For  a time-independent  Hamiltonian,  using  exp(— 2 kH)  = exp(— KH).exp(— kH), 
it  is  easy  to  show  that 


— OO 


(6.6) 
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Starting  with  I\(X,  X';  k),  one  can  iterate  Eq.(6.6)  N times  to  get  I\(X,  X'\  2n k) 
in  terms  of  K(X,X'-,  k).  If  k is  chosen  such  that  k then  N iterations  yield 
the  desired  propagator  in  Eq.(6.5).  If  N is  sufficiently  large  then  k is  small  enough 
to  write  the  short  time  approximation  for  K(X,  X';  k)  as 


-^-(X-X')-I{V(X)  + V(.X')} 

(6.7) 


The  infinite  integral  in  Eq.(6.6)  is  first  replaced  by  a finite  integral;  this  can  be  done 
firstly,  because  the  diagonal  terms  in  Eq.(6.7),  for  which  X'—  X,  are  exponentially 
damped  by  the  contributions  from  the  potential  energy  factor  exp{— kV(.Y)}  for 
large  values  of  X under  the  presumption  that  one  is  dealing  with  binding  potentials, 
and  secondly,  because  the  far  off-diagonal  terms  are  damped  by  the  kinetic  energy 
factor  cxp{—(mx/2h2K)(X—X')2}.  A grid  is  constructed  of  M +1  equally  spaced 
points  over  this  finite  range.  If  A is  the  spacing  between  these  equally  spaced  grid 
points,  and  using  the  trapezoidal  rule  for  the  finite  integral  in  the  coordinate  space, 
each  integration  in  Eq.(6.6)  then  becomes  a simple  numerical  matrix  multiplication 
expressed  as 


M+ 1 

K (i  A,  j A;  2k)  = A ^ K (i A,  /A;  k)  K (/A,  j A;  k)  (6.8) 

/=i 

where  the  real  positive  nature  of  k ensures  the  numerical  stability  of  the  procedure. 

The  advantages  of  using  the  NMM  for  the  integral  in  Eq.(6.6)  are  its  ability 
to  deal  with  any  binding  potentials,  and  the  fast  convergence  towards  the  final 
imaginary  time  ffh  because  during  each  iteration  the  initial  imaginary  time  k 
increases  by  a factor  of  2. 
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The  main  disadvantage  of  the  NMM  method  is  that  the  number  of  grid  points 
needed  to  span  the  appropriate  region  of  a one  dimensional  potential  is  typically  of 
the  order  of  102.  In  one  dimension,  each  NMM  step  for  a 100  x 100  matrix 
requiring  about  106  mathematical  operations,  is  very  well  within  the  limits  of 
present  day  computers.  However,  if  one  were  simply  to  extend  this  procedure 
to  a two  dimensional  problem,  the  matrix  size  of  104  x 104,  with  each  NMM  step 
requiring  roughly  1012  operations,  starts  to  reach  the  limits  of  present  day  computing 
capabilities.  Therefore,  a straightforward  extension  of  the  NMM  method  to  physical 
problems,  involving  more  than  one  or  two  degrees  of  freedom,  is  impractical.  One 
suggested  way  of  improving  upon  this  is  to  use  some  higher  order  quadrature 
instead  of  the  simple  trapezoidal  rule.  Given  the  appropriate  weight  factors  for 
the  selected  grid  points  in  the  coordinate  space,  one  uses  a far  smaller  number  of 
the  total  grid  points  in  each  NMM  step.  This  improvement,  however,  has  never 
been  implemented  in  any  numerical  computations  till  now,  and  its  usefulness  is 
more  in  doubt  if  one  is  using  the  NMM  method  for  computing  the  propagator  for 
a time-dependent  Hamiltonian. 

Thirumalai  and  Berne  recommend  the  use  of  the  NMM  method  only  for 
attractive  or  binding  potentials,  because  otherwise  the  number  of  grid  points  to 
effectively  sample  a repulsive  potential  in  the  coordinate  space  becomes  too  large. 
So  that  even  in  one  dimension  this  increase  in  the  number  of  grid  points  in  the 
NMM  procedure  is  reflected  in  a quadratic  decrease  in  the  computational  efficacy 
of  the  overall  method. 

Lastly,  the  fast  convergence  to  the  final  imaginary  time  /?  for  the  problem  is 
achieved  only  for  time-independent  Hamiltonians;  the  method  as  described  above 
does  not  work  for  a time-dependent  Hamiltonian. 
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6.2.  The  Complex  Time  Propagator  for  the  Primary  Region  Dynamics 


The  NMM  method  for  a one  dimensional  time-independent  Hamiltonian,  can 
be  easily  extended  to  computations  in  the  complex  time  plane  by  replacing  in  the 
above  description  the  imaginary  time  — i/3h  by  a complex  time  r — t — iuh  with 
u as  a constant  real  parameter.  Therefore,  the  short  imaginary  time  k is 

replaced  by  the  short  complex  time  e =jp.  The  complex  time  propagator  matrix 
for  one  dimensional  time-independent  Hamiltonian  is  then  written  as 

K(r)  = Z\N-1.K(e)N_1  (6.9) 

where  if  N is  sufficiently  large  and  e sufficiendy  small,  the  short  time  approximation 
for  the  propagator;  A',y(c)  = K(iA,jA,  e),  is  written  as 


Kii(e)=(.£kyD(i-j)exp 


V 


.2 % 
1 2 ht 


(*-i)2  + f (Yi  + Vj) 


(6.10) 


where  A is  the  spacing  between  any  two  successive  grid  points  in  a (M  + 1)  x 
(M  + 1)  matrix  over  a range  R given  by  R = MA,  and  V,  and  V,  are  the  potentials 
at  the  ith  and  the  jth  grid  points  respectively.  The  contributions  from  the  decay 
factor  D(i  — j),  given  by 


D (i  — j)  = exp 


[ 2N~1mxA 2 


(6.11) 


are  equal  to  one  for  the  on-diagonal  terms  for  which  i = j,  and  decay  rapidly  as 
exp(— C(i  — j)2)  for  i / j,  where  the  decay  constant  C depends  upon  the  choice 
of  the  parameters  M and  N.  The  contributions  from  the  off-diagonal  terms  in  the 
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short  time  complex  matrix  K(e)  can  be  increased  or  decreased  by  decreasing  or 
increasing  the  value  of  the  decay  constant  C. 

In  the  TDSCF  approach  the  effective  frequency  operators,  which  govern  the 
time-evolution  in  the  primary  region  dynamics,  are  time-dependent.  Therefore, 
the  NMM  method  in  which  the  initial  complex  time  e increases  by  a factor  of 
2 during  each  iteration,  does  not  work.  However,  the  complex  time-evolution  in 
the  primary  region  amplitude  £(X,  r)  can  still  be  written  in  a discretized  path- 
integral  formulation  by  discretizing  a finite  region  in  the  coordinate  space  within 
which  the  primary  region  amplitude  f (X,  r)  evolves  on  the  upper  potential  energy 
surface.  For  a primary  region  dynamics,  evolving  with  a time-dependent  effective 
frequency  operator,  the  discretized  PI  formulation  is  then  written  as  a time-ordered 
matrix  product 

K(r)  = z\(N-1).K(N-1)(e).K(N-2)(e)...K(0)(e)  (6.12) 

where  e =jy,  and  K^(e)  is  a complex  time  matrix  at  the  1th  step  in  time  with  / 
going  from  0 to  N — 1.  The  progression  from  the  Ith  step  to  the  (/  4-  l)th  step  is 
simply  a finite  complex  element  matrix  multiplication. 

We  note  that  the  complex  time  TCF  Fb  b(r)  and  its  Fourier  transforms 
Fb  b(uj)  are  computed  on  a contour  along  which  the  the  parameter  u,  which 
defines  the  imaginary  component  of  the  complex  time  r = t — iuh,  remains  a 
constant.  As  shown  in  Fig.l,  the  complex  time  TCF  and  its  Fourier  transform  are 
computed  along  a path  C along  which  the  parameter  u remains  a constant  and  the 
real  time  t'  varies  from  0 to  t — ► oo.  To  effectively  compute  the  Fourier  transform 
along  the  path  C,  one  needs  complex  time  TCF  for  many  values  of  t'  between  0 
and  t.  The  complex  time  r'=  t'  - iuh  corresponding  to  each  of  these  values  of  t' 
is  therefore  divided  by  an  integer  N and  the  time-ordered  matrix  multiplication  for 
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the  primary  region  dynamics  is  performed  at  each  instant  r/along  the  path  C' . To 
complete  the  discussion  about  the  computations  of  a complex  time  propagator  for 
the  primary  region  dynamics,  it  is  natural  that  we  now  consider  the  evaluations  of 
the  complex  time  effective  frequency  operators  along  the  path  C'. 

The  time-dependence  in  the  effective  frequency  operator  Qp(r')  along  a path 
C'  in  the  complex  time  plane,  as  described  in  equations  (5.3)  and  (5.15)  of  the 
previous  chapter,  is  through  the  coupling  between  the  primary  and  secondary 
region  dynamics.  For  simplicity,  we  consider  intermolecular  processes  in  which 
the  effective  coupling  at  the  Ith  step  along  C',  is  written  as 

9 (Tl)  = 9 (0) cos  (WTz)  + ^ sin  (wr()  + y (r/)  (6.13a) 

mu 

where  the  memory  term  y(r/)  is  an  integral  in  the  complex  time  plane 

Tl 

V (Tl)  = [ dtc  v (te)  shuu  (r/  - tc)  (6.13 b) 

mu  J 
o 

and  at  the  Ith  step  r/  = le'  with  e'  = jj.  Provided  the  time-dependent  effective  force 
v(tc)  is  an  analytic  function  in  the  complex  time  plane,  the  value  of  the  complex 
integral  in  the  memory  term  is  independent  of  the  path  C'.  Among  many  possible 
approaches  to  the  evaluation  of  the  integral  in  the  complex  time  plane  we  have 
successfully  used  two  alternate  ways  to  compute  the  memory  term  in  Eq.(6.13b). 
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Figure  1.  Computations  of  the  primary  region  propagator  in  the  complex  time 
plane. 
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The  first  approach  uses  the  analytic  property  of  the  memory  term  in  the  complex 
time  plane;  according  to  which,  the  memory  term  is  first  computed  along  the 
negative  imaginary  time  tc  = — iu'h  for  many  different  values  of  0 < u'  < u, 
and  then  uses  these  values  at  u'  = m to  analytically  continuing  the  memory 
term  to  the  complex  time  plane  at  t\  = t\  — iufi  by  using  a numerical  Pade’ 
approximant  (Schlessinger’s  Point  Method  [Appendix  B])  scheme.  Numerical 
analytic  continuation  for  the  complex  time  computations  has  also  been  used  for 
computing  the  Flux-Flux  TCFs  [Miller  et  al.,  83]  for  reaction  rates  , and  the  dipole- 
dipole  TCFs  [Thirumalai  and  Berne;  83]  for  spectroscopic  properties.  However, 
the  numerical  analytic  continuation  methods  though  accurate  for  short  real  time 
processes,  with  t < uh,  are  not  very  accurate  for  physical  processes  involving  long 
real  times  t » uh.  Therefore,  in  our  second  approach  to  the  computations  of  the 
integral  in  Eq.(6.13b),  we  avoid  the  numerical  analytic  continuation  and  write  the 
integration  of  the  memory  term  in  the  complex  time  plane  as  a sum  of  the  integrals 
along  the  paths  C/  and  the  path  Cf,  i.e 


= J (...)  dtc  + j (...)  dtc 
Cl  Cf 


(6.14) 


where  as  shown  in  Fig.l,  the  path  C\  at  the  real  time  t = 0,  along  the  negative 
imaginary  axis  involves  a real  integral  from  0 — ► —iu\h,  and  is  followed  by  the 
path  Cj2  along  which  the  imaginary  time  —info  remains  constant  while  the  real 
time  varies  from  0 — > </.  Substituting  in  the  first  term  tc  = —iu'h,  and  in  the 
second  term  tc  = t1  - iuih,  the  memory  term  y(r/)  in  Eq.(6.13b)  at  the  complex 
time  77  = ti  — iu\h  is  written  as 
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y(Tl)  — [ du'  v [— iu'h ) sinhu;  (tq  — u')  h 

mu  J 
o 

U 

[ dt'  v [t'  — iu\h)  sinu;  [ti  — t ')  (6.15) 

mw  J 


0 

and  the  integrals  in  the  above  equation  are  evaluated  numerically.  Using  the  same 
procedure  in  the  complex  integrals  in  Eq.(5.20),  we  can  also  simplify  the  propagator 
for  the  secondary  region  dynamics  as 


Ks  [Q,  Q',  t')  = ( — * exp[ 
v ’ \2nih  sin  ut'  J 1 


imu 


2h  sin  u;r' 


x { (q2  + Q'7)  cos  lot'  - 2QQ'  - — /i  (r') 
V / mu  ' 


- ^-h  (r‘) 


*3  (r')  } 


+ { (m  - EY)  r'  ] 


(6.16a) 


where 


U 

I\  [t')  = — J du  v [—iu'h)  sinh  ( uu'h ) 


i 

+ J dt"  v (t"  — iuh)  sin  (t"  — uuh)  (6.166) 


U 

I'2  ( t ')  — — J du'  v [—iu'h)  sinhu;  [u  — u')  h 


i 

+ J dt"  v [t"  — iuh)  sinu;  [t  — t")  (6. 


16c) 
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and 


U u 

J3  (V)  = J duv  (— iu'h ) sinhu;  (u  — u ) h j du"v  (— iu"h ) sinhwu"^ 


o 

U 


— J dt"v  (t"  — iuh)  sinu;  (i'  — t")  J du'v  (— iu'h ) sinhuu/^ 
o o 

t' 

+ J dt"v  ( t " — iuh ) sinu;  (t1  — t") 
o 

t" 

x J dt"'v  (t'"  — iu'h ) sinu;  (tfw  — (6.16d) 


The  complex  time  computations  of  the  primary  region  dynamics  for  a dipole-dipole 
TCF  and  its  Fourier  transform  are  illustrated  in  Fig.  2.  The  dipole-dipole  TCF 
Fjj  £>(t')  is  computed  along  C for  many  points  t'  = t'  — iuh.  The  path-integral 
method,  through  an  efficient  restricted  matrix  multiplication  (RMM)  method  (to  be 
described  in  the  next  paragraph)  for  the  primary  region  dynamics,  is  used  along  the 
path  C'  joining  0 and  r'.  The  memory  terms  needed  to  perform  the  RMM  procedure 
along  C'  are  computed  numerically  according  to  Eq.(6.15)  by  using  the  force  term 
v(t ) at  the  grid  points  shown  in  Fig.2.  The  grid  points  needed  to  compute  the 
memory  term  at  any  complex  time  r/  are  the  points  for  which  the  corresponding 
complex  times  are  less  than  r/;  therefore,  the  force  terms  at  these  grid  points  are 
already  known  from  the  previous  computations. 
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Figure  2.  Complex  time  computations  for  the  TCF  and  its  Fourier  transform; 
discretization  of  the  complex  time  plane. 
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Lastly,  we  explain  the  efficient  Restricted  Matrix  Multiplication  (RMM)  method 
mentioned  in  the  preceding  paragraph.  As  mentioned  before,  even  in  a one 
dimensional  problem,  the  straightforward  numerical  matrix  multiplication  (NMM) 
[Thirumalai  and  Berne,  83]  is  suitable  mostly  for  binding  or  attractive  potentials 
where  the  range  of  the  matrix  sampling  the  potential  is  fixed  by  the  range  of 
the  potential.  Therefore,  for  the  nonbinding  and  repulsive  potentials,  the  method 
becomes  unsuitable  even  in  one  dimensional  problems.  We  use  two  properties  of 
the  short  time  approximation  to  the  complex  matrix  K(e)  in  Eq.(6.10)  to  efficiently 
incorporate  a linear  increase  with  the  range  of  the  potential. 

Firstly,  the  short  time  approximation  for  the  complex  time  matrix  element 
Kij(e)  as  written  in  Eq.(6.10)  is  invariant  to  the  interchange  of  i and  j,  i.e 
Kij(e)  = Kji(e).  Therefore,  by  using  the  property  that  the  product  of  two 
symmetric  matrices  is  a symmetric  matrix,  we  reduce  the  comptational  efforts  in 
any  NMM  iteration  by  a factor  of  two. 

More  importantly,  the  second  property  comes  from  separating,  in  the  short 
time  approximation  for  the  complex  time  matrix  element  K,j(e),  a completly  real 
decay  factor  D(i  — j).  For  diagonal  elements,  i.e  for  j = i,  the  decay  factor 
D(0)  is  one.  For  j i the  decay  factor,  D(i  - j ) oc  exp(-C(i  - j)2),  starts  to 
cause  a rapid  decay  in  the  magnitudes  of  the  off-diagonal  matrix  elements  as  one 
goes  away  from  the  diagonal.  Since  suitable  grid  spacings,  are  generally  achieved 
by  keeping  the  decay  constant  C within  a range  0.1  < C < 1.0  [Thirumalai, 
Bruskin  and  Berne,  83;  Coalson,  85],  the  corresponding  decay  factor  for  far  off- 
diagonal  matrix  elements,  i.e  (i  — j)  > 15  is  typically  of  the  order  of  ~ 10~40. 
This  means  that  the  contributions  from  the  far  off-diagonal  matrix  elements  rapidly 
approach  zero  as  one  goes  further  away  from  the  diagonal.  Instead  of  using  the 
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full  matrix  in  an  NMM  step,  we  use  the  rapidly  decaying  property  of  the  decay 
factor  D(j  — i ),  and  suggest  a natural  cut-off  in  the  contributions  from  the  far 
off-diagonal  matrix  elements.  The  name  Restricted  Matrix  Multiplication  (RMM) 
hence  refer  to  the  NMM  procedure  where  the  contributions  from  the  far  off-diagonal 
elements  beyond  a suitable  cut-off  are  chosen  as  zero.  If  Mc  is  the  suitable  cut 
off  beyond  which  we  are  ignoring  the  contributions  from  the  off-diagonal  elements 
in  a (M  + 1)  x (M  + 1)  matrix,  the  number  of  mathematical  operations  needed 
to  perform  a restricted  matrix  multiplication  (RMM)  step  are  typically  of  the  order 
of  ~ (M  + 1)  x M3  as  opposed  to  the  operations  ~ (M  + l)3  needed  in  an 
equivalent  NMM  step.  Therefore,  as  the  size  of  the  matrix  starts  to  increase  from 
(A/+l)x(M+l)  to  (A/+-L+l)x(M+£+l),  the  RMM  step  for  the  computations 
of  the  same  accuracy  needs  only  ~ L x M3  extra  operations,  whereas  for  the  similar 
increase  in  the  matrix  size  the  NMM  step  will  need  L3  extra  operations.  In  typical 
applications  both  M and  L are  generally  of  the  order  of  ~ 100  whereas  the  cut-off 
parameter  Mc  is  generally  of  the  order  of  ~ 10.  Therefore,  specially  for  unbound 
potentials,  the  efficacy  of  the  RMM  procedure  is  manyfold  as  compared  to  that  of 
the  NMM  procedure. 

6.3.  A More  Efficient  Transformation  from  Real  to  Complex  Time 


As  mentioned  before  the  transformation  from  real  to  complex  time  along  a con- 
tour for  which  the  imaginary  part  -iuh  remains  a constant,  is  not  unique.  We  also 
note  that  having  separate  computation  paths,  for  the  complex  time  Fourier  trans- 
form along  C and  for  the  RMM  steps  along  C',  is  not  very  efficient.  This  is  because 
the  computations  being  done  along  two  separate  contours  require  information  on 
many  more  points  on  the  complex  time  grid,  as  compared  to  a situation  where  C 
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and  C'  coincide,  i.e  where  the  Fourier  transform  computation  is  done  along  the 
RMM  computation  path.  By  deforming  the  path  C for  the  complex  time  Fourier 
transform,  we  show  that  there  exists  another  transformation,  for  which  the  Fourier 
transform  contour  overlaps  with  the  RMM  contour,  making  the  computations  for 
this  new  transformation  more  efficient. 

The  line  shape  function  for  a general  two  surface  photointeraction  process  is 
related  to  the  complex  time  Fourier  transform 


— OO 


(6.17) 


of  the  complex  time  dipole-dipole  TCF  F(t')  with  r'  = t'  — iuh.  We  can  rewrite 


Eq.(6.17)  as 


Fbb  (w)  = lim  dt'  eW-iuh)Fbb  (f  - iuh ) (6.18) 

o 

t — ► OO 


where  we  have  used  F(-t'  - iuh ) =F(t'  - iuh)*,  and  for  constant  u we  have 
taken  the  exponential  factor  containing  u inside  the  integral.  As  shown  in  Fig.(3), 
deforming  the  integration  path  C we  rewrite  the  integral  as 


/ (*')  = / ^‘'hb 

c c, 


+ I dz'e^’Fbt>  (z') 
Cc 


(6.19) 


where  the  complex  integration  variable  z1  = t'  — iuh. 
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Figure  3.  Deformation  of  the  Fourier  transform  contour. 
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In  the  first  term  we  use  z'  = — iu'h  to  show  that 

U 

J dz'e'“*'FDb  {z')  = ih  J du'  e"u'Fbi)  {-iu'h)  (6.20) 

Ci  0 

where  Fb  b(-iu'h)  is  real;  therefore,  the  contribution  to  the  line  shape  function 
from  the  first  term  in  the  deformed  path,  C\ , is  equal  to  zero.  In  the  second  term, 
we  use  r - Vt2  + u2  and  a =arcsin  ■ with  r ->  oo  for  constant  small  angle 
a,  corresponding  to  t -*  oo  for  constant  u,  and  we  rewrite  the  second  term  as 

re“'Q 

J dz'eiuI' F ftp  (z1)  = J dz'  (6.21) 

Cc  0 

Replacing  the  integration  variable  z'  = r'  exp(— ia)  for  constant  a,  and  rearranging 
the  integrand  we  get 

r 

Fib(u)=  lim  ^ e-“  J dr'  (r'e~ia)  (6.22) 

0 

r -t  oo 

where  the  Fourier  transform  integration  and  the  complex  time  dipole-dipole  TCF 
F(r'  exp(— ia))  for  a constant  positive  angle  a,  are  evaluated  along  the  same  path 
Cc.  We  also  note  that,  instead  of  first  translating  the  integration  range  in  the  lower 
half  of  the  complex  time  plane  by  a constant  amount  -iuh,  and  then  deforming 
the  resultant  path  as  C = C\+Cc  to  achieve  the  new  integration  path  in  Eq.(6.22), 
we  could  have  directly  obtained  the  desired  transformation  of  the  integrals  merely 
by  rotating  the  real  time  axis  by  a constant  negative  angle  a = - arc  sin  jt2u+u2 
where  u is  a real  parameter. 


CHAPTER  7 

APPLICATIONS  AND  RESULTS 


In  this  chapter  we  present  applications  of  the  TCF  approach  to  the  interaction 
of  light  with  a polyatomic  system.  Section  7.1  describes  application  to  the  elec- 
tronic absorption  spectra  for  a model  one-dimensional  problem  with  two  electronic 
potential  surfaces.  The  upper  and  the  lower  potential  energy  surfaces  are  modeled 
by  displaced  Morse  potentials  [Coalson,  85]  in  one  dimension.  Computations  of 
the  complex  time  TCF  and  its  Fourier  transform  for  this  model  are  performed  to 
provide  a check  of  the  convergence  of  the  numerical  path-integral  method.  In  sec- 
tion 7.2,  we  use  the  dipole-dipole  TCF  of  an  extended  molecular  system  within  a 
TDSCF  approximation  to  study  the  photodissociation  dynamics  of  methyl  iodide 
(CH3I)  on  the  potential  energy  surface  of  an  electronically  excited  state.  Realistic 
potentials  energy  surfaces  for  this  well  known  example  are  taken  from  a Gaus- 
sian Wave-packet  Dynamics  (GWD)  study  of  the  same  problem  [Lee  and  Heller, 
82].  Using  the  TDSCF  approximation  in  chapters  4 and  5,  the  photodissociation 
dynamics  of  CH3I  has  been  factored  into  a primary  and  a secondary  region  dy- 
namics. The  complex-time  evolution  of  the  primary  region  amplitude  is  performed 
through  an  efficient  discretized  path-integral  method,  whereas  the  complex  time 
propagator  for  the  secondary  region  amplitude  is  constructed  analytically  accord- 
ing to  Eqs.(6.16a,b,c,d)  in  the  previous  chapter.  The  total  cross-section  for  the 
photodissociation  of  CH3I  from  its  ground  vibrational  state,  using  our  TDSCF  ap- 
proach, is  compared  to  that  of  the  GWD  approach  of  Lee  and  Heller.  The  TDSCF 
approach  to  photointeraction  processes  in  many  dimension,  coupled  with  a path- 
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integral  method  for  the  primary  region  dynamics,  is  independent  of  any  particular 
functional  form  of  the  initial  amplitudes  for  the  primary  and  the  secondary  region 
dynamics,  and  of  any  particular  functional  form  for  the  primary  region  potential. 
An  additional  application  to  the  photodissociation  dynamics  of  vibrationally  ex- 
cited CH3I,  for  which  the  initial  values  of  the  primary  and  the  secondary  region 
amplitudes  are  not  restricted  to  simple  Gaussian  functions,  are  described  for  many 
vibrationally  excited  states  in  section  7.3. 

7.1.  Electronic  Absorption  Spectra  of  a Model  Two-Potential  Problem 

For  a two  surface  electronic  absorption  problem  the  complex  time  dipole-dipole 
TCF  t ) of  Eq.(3.11b)  at  the  temperature  (^.g/?)-1  is  written  in  the  coordinate 
representation  as 


energy  surfaces,  the  complex  time  r = t — iuh,  and  Zp  is  the  partition  function 
at  the  initial  inverse  temperature  kgf3.  Electronic  excitations  in  a single  nuclear 
variable  problem,  such  as  in  diatomics,  are  described  by  one  dimensional  potential 
energy  curves.  We  model  the  lower  and  the  upper  potential  energy  curves  as  two 
displaced  one  dimensional  Morse  oscillator  potentials  [Reimers,  Wilson,  and  Heller, 
83;  Thirumalai  and  Berne,  85;  Coalson;  85],  expressed  as 


D (X')  (x'  e~*A,r  D (X)  (7.1) 


where  Hi  and  Hu  are  the  Hamiltonians  for  the  lower  and  the  upper  potential 


(7.2a) 
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Upper  and  Lower  Potential  Energy  Surfaces 


Figure  4.  Model  one  dimensional  Morse  potentials  for  the  ground  , and 

the  excited electronic  potentials  for  a two  state  electronic  absorption 

problem  (see  text). 
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2 

Vu  (X)  = Vo  (l  - e°(A-a))“  - Vo  (1  - eoa)~  . (7.2 b) 

In  what  follows  we  choose  the  same  parameters  as  Coalson’s  [Coalson,  85],  which 
are  listed  in  table  7-1. 

Table  7-1 

Input  Parameters  For  Electronic  Absorption  Problem 

a = 3.0au 
Vo  = 100.25att 
a = 0.071au 
m = l.Oau 
/?  = 0.64au 
u = 0.32au 


The  upper  and  the  lower  potential  energy  surfaces  of  Eq.(7.2a,b)  when  plotted 
in  the  coordinate  space  are  given  in  Fig.4.  Using  (3  = 0.64au  and  the  parameter 
u = 0.32 au  in  Eq.(7.1),  we  rewrite  it  as 


fbbW  = -kJdxJdx'(x 


-Hit- 


X' 


D(X')^X'\e-6'‘T^XSjD(X)  (7.3) 


where  r = t — z0.32,  and  we  have  also  dropped  the  factor  ^ from  the  TCF  because 
it  can  be  taken  care  of  in  the  normalization  of  the  resulting  cross-section.  The 
complex  matrix  for  the  path-integral  method  with  time-independent  Hamiltonian  is 
set  up  by  discretizing  the  coordinate  X within  in  a range  —11  > X > 14  on  a 
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(M+l)x(M+l)  matrix,  and  discretizing  the  complex  time  r = t — i0.32  into  short 
time  steps  of  width  e = The  real  time  t is  varied  between  0 > t > 10  in  the 
incrcaments  of  At  = O.lau.  We  compute  the  upper  surface  propagator  for  many 
values  of  r = t — i0.32  by  N NMM  steps  for  N = 1 with  M = 100,  and  for  AT  = 4 
with  M = 125.  These  choices  of  the  parameters  M and  N keep  the  decay  constant 
C,  as  described  in  chapter  6,  within  the  suitable  range  0.1  < C < 1.0.  As  shown 
in  Eq.(7.3a,b),  the  lower  potential  energy  surface  has  the  same  analytic  form  as  the 
upper  potential  energy  surface,  but  is  displaced  by  a constant  amount  a with  respect 
to  the  upper  potential  energy  suface.  The  lower  surface  propagator  is  computed 
from  the  upper  surface  propagator  by  complex  conjugating  and  displacing  the  matrix 
elements  of  the  upper  surface  propagator  by  an  integer  amount  J = ■£,  where  a is 
the  displacement  between  the  upper  and  the  lower  potential  energy  surfaces,  and 
A is  the  spacing  between  any  two  successive  grid  points.  The  propagators  for  the 
upper  and  the  lower  surface  dynamics  are  then  used  in  Eq.(7.3)  to  compute  the 
complex  time  dipole-dipole  TCF  f ^ ^(t),  the  real  and  the  imaginary  part  of  which 
are  plotted  in  figures  5a  and  5b,  respectively,  as  a function  of  real  time  t. 

Using  the  relation  ^(-t)  = the  real  and  the  imaginary  parts  of 

the  TCF  ^(t),  as  a function  of  positive  real  time  t > 0,  are  sufficient  for  the 
computation  of  the  Fourier  transform  f ^ ^(u)  and  the  corresponding  line  shape 
function,  L( u)  = exp(uu>)f ^ which  is  shown  in  figures  6 and  7.  The 

structure  less  dashed  line  in  Fig. 6,  refers  to  the  Fourier  transform  integral  evaluated 
only  for  short  time  processes,  i.e  the  real  time  limits  in  the  Fourier  transform  integral 
such  as  to  truncate  the  TCF  //>/>(<)  after  its  initial  decay  at  t = 2.0 au. 
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Real  Part  of  the  TCF  vs  Time 


Figure  5a.  Real  part  of  the  complex  time  TCF  //s/s(f)  as  a function  of  real 
time  t; for  N = 4,  M = 125  and for  N = 1,  M = 100. 


Imaginary  Part  of  the  TCF 
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Imaginary  Part  of  the  TCF  vs  Time 


Figure  5b.  Imaginary  part  of  the  complex  time  TCF  /a  t ) as  function  of 

real  time  t; for  N = 4,  M — 125,  and for  N = 1,  M = 100. 


Line  Shape  Function 
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Line  Shape  Function  vs  Frequency 


Figure  6.  Line  shape  function  for  a two-potential  electronic  absorption  problem; 

corresponds  to  retaining  the  short  time  behaviour  in  the  TCF, 

corresponds  to  intermediate  time  behaviour  in  the  TCF  (See  text) 
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Line  Shape  L(W)  vs  Frequency  W 


Figure  7.  Higher  order  overtones  in  the  lineshape  functions  corresponding  to 
long  time  recurrences  in  the  TCF 
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The  highly  oscillating  solid  line  in  Fig.6  corresponds  to  retaining  the  first  recurrence 
in  the  intermediate  time  behaviour,  while  suppressing  any  further  recurrences  in 
the  long  time  behaviour.  This  is  achieved  simply  by  introducing  in  the  Fourier 
transform  integral  a Gaussian  decay  factor,  Ld(t)  — exp (—dt2),  where  d = 0.02. 
If  the  decay  factor  Ld(t)  is  removed  from  the  Fourier  transform  integral,  the  higher 
order  recurrences  in  the  TCF,  which  are  not  even  seen  in  figures  5a  and  5b,  are 
reflected  in  the  lineshape  function  as  the  higher  order  overtone  structures  shown 
in  Fig.7. 

After  checking  the  convergence  of  our  results  by  using  the  full  NMM  procedure 
with  N = 5,  which  agreed  with  the  results  for  N = 4,  we  also  found  agreement 
with  the  results  of  Coalson,  who  has  used  the  same  parameters  and  the  same  method. 
We  repeated  the  computations  for  several  fixed  times,  t = 1.0, 2.0, 3.5,  and  5.0, 
successively  ignoring  more  and  more  far  off-diagonal  matrix  elements  each  time, 
until  a suitable  cut  off  for  the  RMM  method  was  achieved.  We  found  that  for 
N = 4 and  M = 125,  and  Mc  > 25  the  values  of  the  real  and  the  imaginary 
parts  of  the  TCF  using  the  RMM  procedure  remain  converged  to  within  10-6  of 
the  full  NMM  results. 


7.2.  Photodissociation  of  CH3I  from  its  Ground  Vibrational  State. 


We  consider  CH3I  as  a realistic  polyatomic  system  which  undergoes  rearrange- 
ment in  an  electronically  excited  state,  following  the  absorption  of  light  in  the 
visible  and  UV  region.  We  use  a)  the  TDSCF  approximation  for  factoring  the  TCF 
of  the  system  into  self-consistent  primary  and  secondary  region  TCFs,  b)  the  RMM 
iteration  procedure  for  the  primary  region  dynamics,  c)  analytic  construction  of  the 
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secondary  region  dynamics,  and  d)  compute  the  total  TCF  for  the  process,  which 
is  then  used  for  computing  the  total  cross-section. 

The  photodissociation  of  CH3I  molecule  in  an  electronically  excited  state  is 
theoretically  interesting,  firstly,  because  of  the  apparent  validity  of  its  treatment 
as  a pseudotriatomic  molecule  [Shapiro  and  Bersohn,  80;  Lee  and  Heller,  82;  and 
Sundberg  et  al.,  86],  and  secondly,  because  of  its  utilization  as  a test  case  in  many 
of  the  studies  of  photodissociation  in  polyatomic  systems.  In  the  pseudotriatomic 
model,  the  dissociation  of  CH3I  is  described  by  two  coordinates;  rci  is  the  distance 
between  the  C atom  and  the  I atom,  and  vch3  is  the  distance  from  the  C atom 
to  the  plane  containing  the  three  H atoms.  The  nuclear  vibrational  motions  in  the 
molecule  along  the  coordinates  rci  and  rcH3  will  be  refered  to  as  the  Cl  streatch, 
and  the  umbrella  mode  oscillations  of  CH3  respectively.  The  kinetic  energy,  and 
the  ground  and  the  first  excited  potential  energy  surfaces  in  these  two  coordinates, 
are  given  by  Shapiro  and  Bersohn  [Shapiro  and  Bersohn,  80].  An  alternative  set 
of  coordinates,  which  avoid  the  momentum  coupling  in  the  kinetic  energy  terms, 
are  the  mass  weighted  Jacobi  coordinates  [Lee  and  Heller,  82]  R and  r,  given  by 


r = rCH3  (7.4a) 

R = rci  H 7^ rCH3  (7.46) 

The  ground  and  the  first  excited  potential  energy  surfaces  of  CH3I  in  the  Jacobi 
coordinates  are  as  shown  in  Fig.  8.  In  this  work,  however,  we  use  normal  mode 
coordinates  of  CH3I  in  its  ground  electronic  state  [Lee  and  Heller,  82].  The 
transformations  from  the  Jacobi  coordinates  r and  R to  the  mass  weighted  normal 
mode  coordinates  X and  Y are  given  by 
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7.830  -0.176 
0.618  4.939 


)( 


7?  — 4.1670 
r - 0.6197 


(7.5) 


i.e  the  coordinate  X is  dominated  by  the  dissociation  variable  R and  the  coordinate 
F is  dominated  by  the  umbrella  mode  variable  r.  The  normal  mode  analysis  is 
useful  in  assigning  two  oscillator  vibrational  quantum  numbers  (nx,ny)  to  the 
initial  nuclear  vibrational  states  on  the  ground  electronic  potential  energy  surface, 
and  in  providing  harmonic  oscillator  basis  functions  for  the  SCF  approximation 
to  the  initial  nuclear  vibrational  states.  The  ground  and  the  first  excited  potential 
energy  surfaces,  in  the  normal  mode  coordinates  X and  Y,  are  expressed  as 


VgT(X,Y)  = -0.245  + 0.178  [exp  (-0.8995  (X,F))  - 1] 


+ [0.0362  + 0.1101exp(— 0.49147?  ( X , F))] 
x [0.2019F  - 4.5434  x 10"3X 


(7.6a) 


+ 0.6197(1  - exp  (—0.49145  (X,F)))j 


with 


B (X,  Y)  = 0.1287 A'  - 2.4659  x l(r2V 


(7.66) 


and 


Vex  ( X , Y)  = 4.071  x 10_2eip  (-0.1567X  + 0.4327  x 10_1F) 


+ 5.631  x 10_2exp  (-0.1783X  - 0.6044  x 10"2F) 
+ 3.620  x 10_2(— 0.1594  x 10_1JV 


+ 0.2091F  + 0.6197)2  (7.7) 


respectively,  where  all  energies  are  in  Rydberg  and  the  distances  are  in  Bohr. 
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Figure  8.  Two-dimensional  contour  plots  of  the  ground  and  the  first  excited 
electronic  potential  energy  surfaces  of  CH3I  in  Jacobi  coordinates  r and  R. 
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The  ground  state  Hamiltonian,  expanded  to  quadratic  terms  about  the  minimum  in 
the  potential  [Lee  and  Heller,  82],  is  expressed  in  the  normal  mode  coordinates  as 


= 0.4982  x 10-»  (-^  + i*2) 

+ 0.1103  x 10-1  (-lJL  + lr2)+...  , (7.8) 

This  defines  effective  masses,  and  the  force  constants  for  the  harmonic  oscillator 
basis  functions  g,(X)  and  gj(Y ) along  the  normal  mode  coordinates  X and  Y, 
respectively.  The  general  form  of  these  basis  functions  [Messiah,  62]  are  given 
for  example  by 


9i 


mxuxX 2 
2 h 


(7.9) 


where  m\  and  wj,  the  effective  mass  and  the  frequency  for  the  X normal  mode 
respectively,  are  given  by  Eq.(7.8).  Treating  the  normal  mode  coordinate  X as 
the  primary  region  variable,  and  the  normal  coordinate  Y as  the  secondary  region 
variable,  we  expand  the  SCF  solutions  (X  | and  (Y  1 1 pj),  for  the  initial  nuclear 
vibrational  wavefunction  on  the  ground  potential  energy  surface,  as 


(X\tf)  = Y,CiS.(X)  (7.10a) 

» 

w?>  = Ec-?»(y)  - (7-106) 

j 

where  the  expansion  coefficients  C?  and  Cj  are  obtained  by  diagonalizing  the 
full  ground  state  Hamiltonian  with  the  full  initial  nuclear  vibrational  wavefunction 
(XY\ipj)  = (X\ip^)  (Y\jpj).  The  ground  and  first  few  excited  eigenvalues  on 
the  ground  electronic  surface,  obtained  by  using  three  basis  functions  each  in  the 
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expansions  in  Eqs.(7.10a)  and  (7.10b),  are  compared  with  those  of  Lee  and  Heller 
in  the  Table  7-2. 

Table  7-2 

First  Few  Eigenvalues  on  the  Ground  Electronic  Surface 


(nx,ny) 

-Escf 

-E'exact 

(0,0) 

0.109437eV 

0. 108656c V 

(1,0) 

0.179039eL 

0.175237eV 

(2,0) 

0.242616eV 

0.240799eV 

(0,1) 

0.260396eV 

0.260514eV 

(see  APPENDIX  C) 


the  ground  nuclear  vibrational  wave  function  corresponding  to  the  eigenvalue  (0, 0), 
given  by 


V>(0,o)  (X,Y)  = [0.996726<7o  (X)  + 8.07921  x 10"2<7i  (X) 

+ 3.11935  x 1O_302  (X)]  x [0.999921^o  (F) 

- 1.23387  x 10~2gi  (Y)  - 2.27367  x lO"3^  (30] 

(7.11) 

is  used  in  the  self-consistent  treatment  of  the  primary  and  secondary  region  dynam- 
ics on  the  excited  potential  energy  surface.  Following  Lee  and  Heller,  we  assume 
that  the  electronic  transition  dipole  moment  operator,  between  the  ground  and  the 
first  excited  potential  energy  surface,  is  a constant.  Therefore  the  SCF  solutions 
| ipPj)  and  |0]),  on  the  ground  potential  energy  surface  are  equal  to  the  initial  am- 
plitudes £7(0)  and  r/7(0),  respectively.  The  coupling  terms  in  the  excited  state 
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potential  of  Eq.(7.7)  are  expanded  in  powers  of  the  secondary  region  variable  Y, 
and  only  the  terms  which  are  linear  in  Y are  retained.  Therefore,  the  excited  state 
potential  of  Eq.(7.7)  can  be  expressed  in  the  following  simplified  form 

Ves  ( X , Y)  = V°  (X)  + V1  (X)Y  + ±AY2  (7.12a) 

where 


V°  (X)  = 3.620  x 10~2  (0.6197  - 0.1594  x 10_1X)2 

+ 4.071  x l(T2e~0  1567X  + 5.631  x 10~2e_0  1783X 

(7.126) 


V 1 (X)  = 2.0  x 0.2019  x 3.620  x 10~2  (0.6197  - 0.1594  x 10_1A') 
4.071  x io-2e-0  1567*  - 5.631  x 10“2e_0  1783X 


(7.12c) 


and 


A = 2.0  x 3.620  x 10-2  x 0.20192  (7.12d) 

We  have  chosen  to  work  with  this  simplified  form  of  the  excited  state  potential, 
because  we  are  mainly  interested  on  whether  a TDSCF  approximation  can  indeed 
be  used  for  intramolecular  dynamical  processes  involving  dissociation  in  an  elec- 
tronically excited  state,  and  secondly  on,  whether  it  reproduces  the  general  features 
of  the  spectra  of  a dissociating  polyatomic  system.  We  have  found  that  answer  to 


107 


both  these  questions  is  affirmative.  The  analytic  form  of  the  complex  time  prop- 
agator, for  a harmonic  secondary  region  coupled  linearly  to  an  arbitrary  general 
motion  in  the  primary  region  dynamics,  was  described  in  chapter  6. 

To  determine  the  range  of  values  for  the  parameter  u in  the  complex  time 
t = t — iu,  in  atomic  units,  we  first  ignore  the  coupling  between  the  primary  and 
the  secondary  regions,  and  at  t = 0 use  the  full  NMM  method  for  the  primary 
region  dynamics  along  the  negative  imaginary  time  axis  for  different  values  of  u 
in  the  range  0 > it  > 100.  To  do  this,  we  use  imaginary  time  steps  of  length 
e = 3au,  and  a coordinate  grid  spacing  A = 0.106au  on  a 75  x 75  matrix  in  a 
range  -4  > X > +4,  within  which  the  initial  value  of  the  amplitude  is  sufficiently 
greater  than  zero.  These,  together  with  the  effective  mass  mx  = 200.723au  along 
the  normal  mode  coordinate  X,  keep  the  decay  constant,  C = 0.38,  within  the 
range  0.1  < C < 1.0  [Thirumalai  and  Berne,  83;  Coalson,  85].  The  upper  bound 
on  u is  fixed,  firstly,  by  noticing  that  for  u > 40au  the  value  of  the  primary  region 
TCF  at  the  initial  time  t = 0 is  smaller  than  about  10-8,  and  secondly,  by  knowing 
that  for  the  complex  time  dynamics  this  value  will  decrease  even  further.  For  the 
lower  bound  on  u the  complex  time  computations  for  the  primary  region  dynamics 
were  performed  for  u between  1 < u < 15.  We  noticed  that  for  u < 7 the  results 
for  the  complex  time  primary  region  TCF  do  not  remain  bound,  i.e  the  amplitude 
of  the  envelop  of  the  oscillating  TCF  start  growing  to  values  which  are  larger  than 
the  amplitude  of  the  first  maxima.  Consequently,  the  suitable  range  for  the  allowed 
values  of  u is  chosen  as  8 < u < 40.  For  the  set  of  results  presented  in  this  section, 
we  have  used  u = 15au,  and  have  checked  that  these  values  remain  converged  to 
within  10"6  for  varying  u in  the  range  10  < u < 35. 

For  the  complex  time  primary  region  dynamics  (including  the  self-consistent 
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coupling  with  the  secondary  region  dynamics),  we  discretize  the  space  along  the 
normal  mode  coordinate  A'  in  a range  —6  > X > 8,  within  which  the  complex 
time  primary  region  amplitude  £( X , r)  evolves  and  spreads  on  the  upper  potential 
energy  surface,  on  a 115  x 115  matrix.  For  fixed  u = 15au,  the  complex  time 
t = t — iu  is  incremented  from  t = 0 — il5  by  increasing  real  time  t in  steps  of 
At  = 4 au  until  the  total  TCF  is  decayed  sufficiently.  Each  of  these  complex  times 
r is  distretized  into  short  complex  time  steps  e = with  N = 5.  These  choices 
for  the  parameters  M and  N , keep  the  the  decay  constant,  C = 0.505,  within 
the  required  range.  The  complex  time  computations,  using  the  full  NMM  method 
for  several  fixed  complex  times  r were  performed  repeatedly  by  dropping  more 
and  more  far  off-diagonal  elements  from  the  computations,  untill  a suitable  cut-off 
point  was  reached.  For  the  set  of  parameters  chosen  above,  i.e  for  C = 0.505,  a 
cut-off  point  Mc  > 18  keeps  the  results  of  the  RMM  method  converged  to  within 
10-6  of  those  of  the  full  NMM  method. 

Table  7-3 

Computational  Parameters  for  the  Primary  Region  Dynamics 

u = 15au 
A min  = 6a  U 
Xmax  — +8au 
M = 115 
Mc  = 25 
N = 5 

Comparisons  of  the  complex  time  evolution  of  the  primary  and  secondary  region 
amplitudes  are  done  in  figures  9a,b,c  and  10a,b,c;  where  at  different  real  times  t, 
we  plot  the  real  parts  of  the  complex  primary  and  the  secondary  region  amplitudes. 
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Primary  Region  Amplitude 
For  Ground  Vibrational  State  (0,0) 


Figure  9a.  Real  part  of  the  primary  region  amplitude  along  the  normal  mode 

coordinate  X at  the  real  time  t = 0 au , t = 16  au , and  t = 32 

au  
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Primary  Region  Amplitude 
For  Ground  Vibrational  State  (0,0) 


Figure  9b.  Real  part  of  the  primary  region  amplitude  along  the  normal  mode 

coordinate  X at  the  real  time  t = 64  au  ....  , t = 128  au , and  t = 256 

au 
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Primary  Region  Amplitude 
For  Ground  Vibrational  State  (0,0) 


Figure  9c.  Real  part  of  the  primary  region  amplitude  along  the  normal  mode 
coordinate  X at  the  real  time  t = 360  au 
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Secondary  Region  Amplitude 
For  Ground  Vibrational  State  (0,0) 


Figure  10a.  Real  part  of  the  secondary  region  amplitude  along  the  normal 

mode  coordinate  Y at  the  real  time  t = 0 au , t = 16  au , and 

t = 32  au  
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Secondary  Region  Aimjlitude 
For  Ground  Vibrational  State  (0,0) 


-8  -6  -4  —2  y 0 ^ 2 4 6 8 


Figure  10b.  Real  part  of  the  secondary  region  amplitude  along  the  normal 

mode  coordinate  Y at  the  real  time  t = 64  au  ....  , t = 128  au , and 

t = 256  au  . 
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Secondary  region  Amplitude 
For  Ground  Vibrational  State  (0,0) 


Figure  10c.  Real  part  of  the  secondary  region  amplitude  along  the  normal  mode 
coordinate  Y at  the  real  time  t = 360  au 
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For  the  ground  vibrational  (0,0)  state  of  CH3I,  the  amplitudes  for  the  primary 
and  the  secondary  region  dynamics  start  as  simple  gaussian  functions  at  the  ini- 
tial real  time  t = 0 and  the  corresponding  complex  time  r = — iu\  as  the  later 
real  time  t > 0 increases,  these  amplitudes  move  from  their  initial  positions,  and 
also  simultaneously  change  their  shape.  In  particular,  the  contrast  between  the 
fast  oscillating  motion  of  the  primary  region  amplitude  and  the  slowly  oscillating 
motion  of  the  secondary  region  amplitude,  is  clearly  displayed  in  these  figures. 
The  complex  time  amplitudes  r)  and  ^(°>°)(F,  r)  are  then  used  for  com- 

puting the  dimensionless  primary  region  TCF  /(°>°)(r),  the  secondary  region  TCF 
s(°-o)(T),  and  the  time-dependent  phase  factor  a^°’°)(r)  as  defined  in  equations 
(4.29b, c,d).  The  primary  and  the  secondary  region  TCFs  together  with  the  phase 
factor  are  sufficient  to  describe  the  complex  time-dependence  in  the  the  total  TCF, 
F(°’°)(r).  xhe  superscript  (0,0)  in  these  expressions  means  that  on  the  ground 
electronic  surface  there  are  no  excitations  along  the  normal  mode  coordinates  X 
(dominated  by  the  Cl  stretch  mode),  and  none  along  the  normal  mode  coordinate 
Y (dominated  by  the  umbrella  mode  vibrations  in  the  CH3  fragment).  The  real 
parts  of  the  complex  TCFs;  /(0-°)(t),  fif(0’0)(r),  and  F(0’°)(r)  as  function  of  real 
time  t are  plotted  in  figures  11,  12  and  13,  respectively.  Finally,  in  figure  14  we 
compare  the  total  cross-section  for  the  photodissociation  of  CH3I  from  its  ground 
vibrational  state,  plotted  as  a function  of  incident  photon  energy,  with  those  of  the 
Gaussian  Wavepacket  Dynamics  method  [Lee  and  Heller,  82].  The  discrepencies  in 
the  peak  position,  and  in  the  full  width  at  half  maximum  (FWHM)  are  very  small, 
and  may  be  attributed  to  our  approximations  expanding  the  excited  state  potential 
in  the  secondary  region  variable  Y and,  retaining  only  the  terms  linear  in  Y,  and 
using  the  SCF  form  for  the  initial  nuclear  vibrational  state  ^0  0)(X,Y). 
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Real  Part  of  the  Primary  Region  TCF 
For  Ground  Vibrational  State  (0,0) 


Figure  11.  Real  part  of  the  dimensionless  primary  region  TCF  f(0,0\t)  as  a 
function  of  real  time  t. 
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Secondary  Region  TCF 
For  Ground  Vibrational  State  (0,0) 


Figure  12.  Real  pan  of  the  dimensionless  secondary  region  TCF  ff(0,0)(t)  as 
a function  of  real  time  t. 
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Real  Part  of  the  Total  TCF 
For  Ground  Vibrational  State  (0,0) 


Figure  13.  Real  part  of  the  total  TCF  _F(°’°)(t)  as  a function  of  real  time  t. 
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Peak  Normalized  Total  Cross-section 


Figure  14.  Peak  normalized  total  cross  section  for  the  photodissociation  of 

CH3I  from  its  ground  vibrational  state; TDSCF  approach, 

Gaussian  Wavepacket  Dynamics  approach  [Lee  and  Heller;82], 
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7.3.  Photodissociation  of  CH3I  from  initially  excited  vibrational  states. 

The  main  advantage  of  using  our  TDSCF  approach  in  more  than  one  dimension, 
is  the  absence  of  any  constraint  on  the  allowed  functional  form  of  the  time- 
dependent  amplitudes  for  the  primary  and  secondary  region  dynamics,  and  the 
primary  region  potential.  Consequently,  the  extension  of  the  TDSCF  approach 
to  the  photodisociation  dynamics  of  CH3I  in  initially  excited  vibrational  states, 
where  the  initial  form  of  the  complex  amplitudes  are  not  restricted  only  to  simple 
Gaussian  shapes,  is  straightforward.  In  the  above  description  of  the  self-consistent 
primary  and  secondary  region  dynamics,  we  simply  replace  the  SCF  form  of 
the  ground  vibrational  state  V’(o,o)(^>  Y)  by  the  required  excited  vibrational  state 
V’(nXlny)(-^5^)  where  nx  and  ny  are  the  number  of  vibrational  quanta  along  the 
dissociation  dominated  normal  mode  coordinate  X,  and  the  umbrella  mode  motion 
dominated  normal  mode  coordinate  Y,  respectively.  The  peak  position  and  the 
shape  of  the  total  cross-section  in  the  frequency  spectrum  are  determined  by  the 
primary  region  TCF,  whereas  the  inclusion  of  the  secondary  region  TCF  provide 
a modulation  in  the  peak  shape  of  the  cross-section.  As  we  consider  the  inclusion 
of  excited  vibrational  states  ( nx  = 1,2,3...)  in  the  primary  region  dynamics,  we 
find  that  increasing  number  of  nodes  in  the  initial  amplitude  are  directly  reflected 
in  the  appearance  of  extra  peaks  in  the  total  cross-section.  On  the  other  hand,  a 
similar  inclusion  of  excited  vibrational  states  (ny  = 1,2,3...)  along  the  secondary 
region  dynamics  does  not  directly  show  in  extra  peaks  in  the  cross-section.  The 
computational  parameters  for  the  dynamics  in  the  primary  and  the  secondary  region 
are  same  as  those  of  in  the  Table  7-2.  In  figures  (15a,b,c)  and  (16a,b,c)  we  show 
the  time  developments  for  the  primary  and  the  secondary  region  amplitudes  for 
initially  excited  states  (1,0)  and  (0,1),  respectively. 
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Primary  Region  Amplitude 
For  Initially  Excited  (1,0)  State 


Figure  15a.  Real  part  of  the  primary  region  amplitude  for  initially  excited 

(1,0)  state  of  CH3I  at  different  real  times  t = 0 au , t = 16  au 

, and  t = 32  au  


122 


Primary  Region  Amplitude 
For  Initially  Excited  (1,0)  State 


Figure  15b.  Real  part  of  the  primary  region  amplitude  for  initially  excited 

(1,0)  state  of  CH3I  at  different  real  times  t = 64  au  ....  , t = 128  au 

, and  t = 256  au . 
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Primary  Region  Amplitude 
For  Initially  Excited  (1,0)  State 


Figure  15c.  Real  part  of  the  primary  region  amplitude  for  initially  excited  (1,0) 
state  of  CH3I  at  the  real  time  t — 360  au 
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Secondary  Region  Amplitude 
For  Initially  Excited  (0,1)  State 


Figure  16a.  Real  part  of  the  secondary  region  amplitude  for  initially  excited 
(0, 1)  state  of  CH3I  at  t = 0 au , t = 16  au , and  t — 32  au 
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Secondary  Region  Amplitude 
For  Initially  Excited  (0,1)  State 


Figure  16b.  Real  part  of  the  secondary  region  amplitude  for  initially  excited 

(0, 1)  state  of  CH3I  at  t = 64  au  ....  , t = 128  au , and  t = 256 

au  . 
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Secondary  Region  Amplitude 
For  Initially  Excited  (0,1)  State 


Figure  16c.  Real  part  of  the  secondary  region  amplitude  for  initially  excited 
(0, 1)  state  of  CH3I  at  t = 360  au . 
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The  complex  amplitudes  corresponding  to  the  initially  excited  (1,0)  and  (0,1) 
nuclear  vibrational  states  are  then  used  for  computing  the  total  TCFs  F(1,0)(r)  and 
_FH0,1)(t),  respectively.  The  real  parts  of  these  TCFs  as  functions  of  real  time 
t and  the  resulting  total  cross-sections  as  functions  of  photon  energy  are  shown 
in  Figs.  17a, b and  Fig.  18,  respectively.  Modulations  of  the  fast  oscillations  in  the 
TCF,  corresponding  to  the  (1,0)  vibrationally  excited  CH3I  as  shown  in  Fig.  17a, 
are  reflected  in  the  appearance  of  an  additional  peak  in  the  corresponding  total 
cross-section  in  Fig.  18,  whereas  the  modulations  of  the  slow  decay  of  the  TCF 
in  the  (0,1)  case  as  shown  in  Fig.l6b  appear  merely  as  a shoulder  structure  in 
the  corresponding  total  cross-section.  At  present  there  are  no  quantum  mechanical 
computations  for  the  photodissociation  of  CH3I  from  initially  excited  vibrational 
states.  Flowever,  qualitatively,  the  general  features  of  our  total  cross-sections  for 
(1,0)  and  (0,1)  vibrationally  excited  CH3I  agree  well  with  two  semi-classical 
computations,  i.e  the  semi-classical  wave-packet  dynamics  method  [Sundberg  et  al., 
86]  for  differently  parametrized  potential  energy  surfaces,  and  the  Self-consistent 
Eikonal  Method  (SCEM)  [Swaminathan,  Stodden  and  Micha,  88]  for  the  same 
potential  energy  surfaces.  Both  of  these  methods  also  show  two  peaks  in  the  total 
cross-section  of  the  photodissociation  of  the  (1,0)  initially  excited  CH3I,  and  a 
single  non-symmetric  peak  in  the  (0, 1)  case.  Qualitatively,  the  shapes  of  our  total 
cross-sections  are  closer  to  the  wavepacket  computations  [Sundberg  et  al.,  86],  as 
opposed  to  the  SCEM  computations,  which  uses  the  same  potential  energy  surfaces 
as  in  our  work.  This  is  probably  because  in  the  SCEM  technique,  which  is  based 
on  running  classical  trajectories  for  the  nuclear  motions,  the  low  energy  region  of 
the  spectrum  is  not  very  accurate,  because  at  low  energies  the  classical  trajectories 
do  not  properly  sample  the  required  non-allowed  region  of  the  potential. 
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Real  Part  of  the  Total  TCF 
For  Vibrationally  Excited  State  (1,0) 


Figure  17a.  Real  part  of  the  complex  TCF  F(1,0)(r)  as  a function  of  real  time  t. 
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Real  Part  of  the  total  TCF 
For  Vibrationally  Excited  State  (0,1) 


Figure  17b.  Real  part  of  the  complex  TCF  F(0,1)(r)  as  a function  of  real  time  t. 
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Normalized  Cross-sections 
For  Initially  Excited  (1,0)  and  (0,1)  States 


Figure  18.  Total  cross-section  for  vibrationally  excited  CH3I  for  a single 

excitation  in  the  secondary  region  dynamics (0, 1),  and  a single 

excitation  in  the  primary  region  dynamics (1,0). 


131 


Using  five  basis  functions  each  in  the  SCF  expansions  for  (X\rp^)  in  Eq.(7.10a) 
and  ( Y\ipj ) in  Eq.(7.10b),  we  obtained  nuclear  vibrational  wavefunctions  and  the 
corresponding  eigenvalues  for  the  other  excited  states  as  well  [Appendix  C],  and 
have  used  them  in  further  studies  of  the  photodissociation  of  vibrationally  excited 
CH3I.  In  Fig.l9a,  we  keep  the  vibrational  quantum  number  ny  fixed  to  ny  = 0 and 
check  the  effects  of  increasing  the  number  of  vibrational  excitations  along  the  the 
normal  mode  coordinate  X.  We  note  that  for  each  additional  vibrational  excitation 
in  the  primary  region  dynamics,  there  is  an  additional  extra  peak  in  the  spectrum. 
The  heights  of  these  additional  peaks  are  not  equal,  peaks  in  the  high  energy  regions 
of  the  spectrum  corresponding  to  the  short  time  behaviour  of  the  TCF  are  found 
to  be  higher  than  the  peaks  in  the  low  energy  regions  of  the  cross-section.  As 
we  increase  the  number  of  vibrational  quanta  in  the  initial  state,  we  note  that  the 
overall  center  of  the  spectrum,  which  is  a maximum  for  even  number  of  vibrational 
excitations  such  as  in  (0,0)  and  (2,0)  cases,  and  is  a minimum  for  odd  number  of 
vibrational  excitations  such  as  in  (1,0)  and  (3,0)  cases,  gradually  shifts  towards 
the  low  energy  region  of  the  spectrum.  In  Fig.  19b,  we  keep  ny  fixed  to  ny  = 1, 
and  check  the  effects  of  increasing  the  number  of  initial  vibrational  excitations  in 
the  primary  region  dynamics.  In  addition  to  the  appearance  of  a shoulder  structure 
due  to  the  excitation  in  the  secondary  region  dynamics,  we  also  note  a growth  in 
the  size  of  this  structure  as  we  increase  the  number  of  vibrational  excitations  in 
the  primary  region  dynamics,  so  that  in  (2, 1)  and  (3, 1)  vibrationally  excited  CH3I 
this  structure  becomes  like  an  additional  peak.  We  also  note  that  the  heights  of 
the  peaks  in  (1,1),  (2, 1),  and  (3, 1)  vibrationally  excited  CHI  are  generally  bigger 
than  the  heights  of  the  corresponding  peaks  in  the  (1,0),  (2,0),  and  (3,0)  cases. 
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Normalized  Cross-sections 
For  Vibrational  Excitations  In  The  Cl  Mode 


Figure  19a.  Comparison  of  cross  sections  for  vibrational  excitations  in  CH3I 

along  the  Cl  stretch  mode:  (1,0) ,(2,0) , and  (3,0) , with 

the  cross-section  for  no  excitation  in  the  umbrella  mode:  (0,0)  (•). 
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Normalized  Cross-sections 
For  Vibrational  Excitations  In  The  Cl  Mode 


Figure  19b.  Comparison  of  cross  sections  for  vibrational  excitations  in  CH3I 

along  the  Cl  stretch  mode:  (1, 1) , (2, 1) , and  (3, 1) , with 

the  cross-section  for  a single  excitation  in  the  umbrella  mode:  (0, 1)  (•). 
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In  Fig.20a,  we  keep  the  quantum  number  n\  corresponding  to  the  initial  vibrational 
excitations  in  the  primary  region  dynamics  fixed  to  nx  = 0,  and  see  the  effect  of 
increasing  the  number  of  initial  vibrational  excitations  in  the  secondary  region 
dynamics.  As  before,  we  note  the  appearance  of  a shoulder  structure  in  (0,1) 
vibrationally  excited  CH3I,  and  that  this  gradually  grows  almost  to  an  additional 
peak  in  the  (0,3)  case.  The  heights  of  the  main  peak  in  the  (0,1),  (0,2)  and  (0,3) 
peaks  are  also  modulated  as  we  put  more  and  more  initial  vibrational  energy  in  the 
secondary  region  dynamics.  In  Fig.20b,  we  fix  nx  = 1 and  gradually  increase  the 
number  of  vibrational  excitation  energy  in  the  secondary  region  dynamics.  Here 
also  the  two  peak  structure  of  the  (1, 0)  case  is  modulated  by  an  additional  shoulder 
structure  in  the  (1, 1)  case,  which  grows  and  spreads  towards  the  low  energy  region 
of  the  spectrum  in  the  (1,2)  and  (1,3)  cases.  We  also  note  a modulation  in  the 
heights  of  the  main  peaks  as  we  increase  the  number  of  initial  excitations  in  the 
secondary  region  dynamics.  Hence,  in  addition  to  being  able  to  produce  highly 
structured  cross-sections  for  the  photodissociation  of  vibrationally  excited  CH3I, 
as  shown  in  figures  19a,b  and  20a, b;  the  TDSCF  approach  to  the  dynamics  also 
shows  the  effect  of  a self-consistent  energy  redistribution  between  the  primary  and 
the  secondary  regions  in  these  cases. 
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Normalized  Cross-sections 
For  Vibrational  Excitations  In  The  CH3  Mode 


Figure  20a.  Comparison  of  cross  sections  for  vibrational  excitations  in  CH3I 

along  the  umbrella  mode:  (0,1) ,(0,2) , and  (0,3) .with  the 

cross-section  for  zero  excitation  in  the  Cl  stretch  mode:(0,0)  (•). 
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Normalized  Cross-sections 
For  Vibrational  Excitation  In  The  CH3  Mode 


Figure  20b.  Comparison  of  cross  sections  for  vibrational  excitations  in  CH3I 

along  the  umbrella  mode:  (1, 1) , (1,2) , and  (1,3) .with  the 

cross  section  for  a single  excitation  in  the  Cl  stretch  mode  (1,0)  (•). 


CHAPTER  8 

DISCUSSION  AND  CONCLUSIONS 


In  chapters  2 through  6 we  have  developed  a general  time-dependent  quantum 
mechanical  approach  to  the  interaction  of  visible  and  UV  light  with  extended 
molecular  systems.  Light  sources  in  this  frequency  range  generally  excite  nuclear 
vibrational  motions,  accompanied  by  transitions  between  different  electronic  states 
of  the  system.  We  have  considered  the  interaction  of  visible  and  UV  light  with  a 
large  molecular  system  as  a two  step  process:  the  interaction  of  light  with  a large 
molecule  system,  treated  like  a photon-molecule  collision  process,  excites  electronic 
transitions  in  the  molecular  target;  and  is  followed  by  a dynamical  evolution  of  the 
target  on  the  excited  potential  energy  surface.  The  discussion  of  our  work  in  this 
chapter  is  also  accordingly  ordered  , i.e  first,  we  discuss  treatment  of  light-molecule 
interactions  as  a photon  collision  process,  and  then  consider  the  dynamical  evolution 
of  a large  molecular  system,  from  the  given  initial  conditions  to  the  required  final 
conditions.  At  the  end,  we  discuss  application  and  results  of  this  approach  to 
study  the  electronic  absorption  spectra  of  a model  two  potential  problem,  and  the 
photodissociation  dynamics  of  a realistic  molecular  system  (CH3I)  from  its  ground 
and  excited  nuclear  vibrational  states. 

We  have  described  transition  rates  for  all  first  and  second  order  photointeraction 
processes  as  the  time  integrals  of  the  product  of  the  molecular  TCFs  of  the  target, 
and  the  radiation  field  TCFs  of  the  light  source.  The  method  of  separating  the 
light-molecule  interaction  into  a radition  field  factor  and  a molecular  target  factor 
is  completely  general,  and  can  be  extended  to  a similar  separation  of  the  light- 
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molecule  interaction  to  higher  order.  An  explicit  incorporation  of  the  statistical  and 
the  dynamical  properties  of  the  light  source  in  the  formalism  is  achieved  by  the 
separation  of  the  radiation  field  TCFs  from  the  molecular  TCFs.  As  examples  of 
the  inclusion  of  the  statistical  properties  of  the  light  source  in  the  formalism,  we 
have  constructed  the  radiation  field  TCFs  of  thermal,  chaotic  and  coherent  light 
sources.  We  have  found,  that  the  expressions  for  our  transition  rates  for  thermal 
light,  in  the  time  domain  description,  are  equivalent  to  the  usual  energy  domain 
descriptions  [Louisell,  73;  Loudon,  83],  except  that  a statistical  distribution  for 
the  number  of  photons  is  also  included  in  the  transition  rates.  In  the  case  of  the 
chaodc  light  source,  which  is  a laser  operating  below  the  threshold,  the  form  of 
the  transition  rates  are  found  to  be  the  same  as  those  in  the  thermal  light  case. 
The  exact  numbers  of  photons  in  the  transition  rates  for  the  thermal  light  are 
replaced  by  a mean  number  of  photons  in  the  transition  rates  for  chaotic  light.  We 
also  note,  that  this  mean  number  of  photons  is  described  by  the  initial  statistical 
distribution  of  the  electric  field  amplitude  of  the  experimental  laser  source  used.  For 
the  coherent  light-molecule  interaction  processes,  we  have  found  additional  terms  in 
the  transition  rates,  due  to  the  coherent  nature  of  the  light  source.  The  magnitudes 
of  these  additional  transition  rates  are  described  by  the  initial  statistical  distribution 
of  the  electric  field  amplitude  and  phase,  given  by  an  experimental  laser  source. 

In  chapters  3 and  4 we  have  discussed  the  evolution  of  a large  system  through  a 
molecular  TCF  approach.  Using  the  Bom-Oppenheimer  approximation  to  separate 
electronic  and  nuclear  motions  in  the  target,  we  have  written  a general  molecular 
TCF  for  a two-surface  electronic  excitation  problem  in  the  electronic  adiabatic 
representation.  A transformation  from  real  to  complex  time  of  these  molecular  TCFs 
is  also  described,  so  that  their  computations  by  a numerical  path-integral  method 
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become  feasible  . The  transformation  to  the  complex  time  of  the  TCFs,  in  which  the 
imaginary  part  of  the  complex  time  remain  fixed,  is  not  very  efficient.  Therefore, 
in  a later  chapter,  we  have  derived  another  transformation,  which  achieves  the  same 
objective,  and  computationally  is  more  efficient  to  implement.  Furthermore,  we  also 
show  that  the  complex  time  dependence  in  these  TCFs  is  described  by  the  dynamical 
evolution  of  an  amplitude  under  the  influence  of  a frequency  operator.  We  have 
factored  the  amplitude  of  the  total  system  into  self-consistent  primary  and  secondary 
region  amplitudes  and  a time  dependent  phase  factor.  A TDSCF  approximation  for 
this  purpose  is  derived  such  that  norm  conservation  in  the  TDSCF  approximation 
and  a correct  phase  factor,  essential  requirements  in  our  work,  are  also  satisfied.  The 
correct  phase  factor  is  required,  because  in  computing  the  TCFs  from  the  complex 
time  amplitudes,  the  time-dependent  phase  factor  does  not  get  canceled.  The  norm 
conserving  form  for  the  factorization  of  the  amplitude  is  essential,  because  the 
effective  frequency  operators,  (f  j j(t)  and  &stff(T)  which  govern  the  dynamics  of 
the  self-consistent  primary  and  secondary  regions,  respectively,  are  non-hermitian 
for  all  complex  time  r > 0.  The  main  advantage  of  a TDSCF  approximation,  and 
the  factorization  of  a multivariable  dynamics  into  primary  and  secondary  region 
dynamics,  is  the  allowance  of  a self-consistent  energy  exchange  between  strongly 
coupled  primary  and  secondary  regions  during  the  evolution. 

Many  intramolecular  photointeraction  processes  in  large  systems  are  modeled 
by  considering  an  arbitrary  general  motion  in  the  primary  region  dynamics,  strongly 
coupled  to  many  harmonic  degrees  of  freedom  in  the  secondary  region  dynamics. 
The  strong  coupling  between  the  primary  and  the  secondary  region  dynamics  is 
modeled  by  considering  the  terms  in  the  coupling  which  are  completely  general  in 
the  variables  of  the  primary  region,  and  have  a linear  and  bilinear  dependence  in 
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the  variables  of  the  secondary  region.  Weak  coupling  limit  of  this  case,  which  is 
more  suitable  for  intermolecular  processes,  is  also  considered  by  ignoring  the  terms 
in  the  coupling  which  are  bilinear  in  the  variables  of  the  secondary  region.  For  both 
of  these  cases,  the  solutions  to  the  equations  of  motion  for  the  secondary  region 
dynamics  are  reduced  to  simple  analytic  forms.  We  have  used  these  solutions  for 
constructing  the  time-dependent  couplings  to  the  primary  region  dynamics,  and  also 
for  constructing  the  propagators  for  the  secondary  region  dynamics.  The  form  of 
the  time-dependent  couplings  are  such  that  the  self-consistent  solution  of  the  full 
problem  is  possible  without  going  through  many  iterations  at  each  step  in  time. 

We  have  discussed  a direct  relationship  between  our  TDSCF  approach  and 
the  Generalized  Langevin  Equation  (GLE)  approach  to  the  dynamics  of  large 
systems.  An  operator  equation  for  the  primary  region  dynamics,  which  contains 
a fluctuation  force  term  and  a memory  term,  is  derived.  We  have  shown  that 
the  fluctuation  force  term  in  this  operator  equation  is  dependent  upon  the  initial 
conditions  of  the  secondary  region  dynamics  only  through  the  expectation  values 
of  the  secondary  region  operators.  The  fluctuation  force  term  is  therefore  given 
a stochastic  interpretation  by  observing  that  the  initial  values  of  these  expectation 
values  are  described  by  a distribution  function  which  depends  upon  the  initial 
temperature  of  the  system.  If  initially,  at  t = 0,  the  secondary  region  is  brought  to 
a thermal  equilibrium  in  the  presence  of  the  primary  region,  then  we  can  consider  a 
displaced  Boltzmann  distribution  function  for  the  initial  values  of  these  expectation 
values,  and  explicitly  derive  a Fluctuation-Dissipation  Relationship  (FDR)  between 
the  fluctuation  force  and  the  memory  kernel  in  our  TDSCF  approach.  The  use  of  the 
displaced  Boltzman  distribution  function,  instead  of  the  usual  Boltzman  distribution 
function,  is  required  to  satisfy  the  central  limit  theorem  for  the  fluctuation  force  at 
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t = 0,  and  also  to  facilitate  a correct  FDR  between  the  ensemble  averaged  TCF  of 
the  fluctuation  force  and  the  memory  kernel. 

We  have  also  discussed  treatment  of  many  degrees  of  freedom  in  the  primary 
region,  coupled  linearly  and  bilinearly  to  many  harmonic  degrees  of  freedom  in  the 
secondary  region.  We  have  shown  that  for  a few  degrees  of  freedom  in  the  primary 
region,  it  is  possible  to  explicitly  treat  the  primary  region  dynamics  without  further 
approximations. 

Using  a discretized  coordinate  space,  we  have  developed  an  efficient  path- 
integral  method  for  the  primary  region  dynamics  that  we  have  referred  to  as  the 
Restricted  Matrix  Multiplication  (RMM)  method.  The  RMM  method  is  based  on  the 
recently  introduced  numerical  matrix  multiplication  (NMM)  method  [Thirumalai, 
Bruskin  and  Beme,  83],  but  avoids  its  problems,  which  include  a quadratic  scaling 
of  the  computational  efforts  with  respect  to  any  increase  in  the  range,  and  the 
dimensionality  of  the  potentials,  and  the  inability  to  deal  with  time-dependent 
potentials.  The  efficient  RMM  method  uses  matrix  multiplication  between  a 
complex  time  matrix  corresponding  to  the  propagator  up  to  a certain  time  and  a 
short  complex  time  matrix  corresponding  to  the  short  time  propagator  for  the  next 
step  in  time.  Since,  at  each  step  in  time,  we  are  using  the  short  time  propagator  for 
the  previous  step,  it  is  easy  to  deal  with  the  time-dependent  potentials  in  the  RMM 
method.  The  term  restricted  matrix  multiplication  refers  to  the  fact  that  instead  of 
using  full  matrices  in  the  multiplication  at  any  given  step  in  time,  we  have  used 
the  property  that  the  off-diagonal  matrix  elements  in  the  complex  time  propagator 
decrease  rapidly  away  from  the  diagonal,  to  find  a suitable  cut-off  point.  The 
contributions  of  the  far  off-diagonal  matrix  elements  beyond  this  cut-off  point,  in 
the  RMM  method,  are  ignored  without  affecting  the  accuracies  of  the  computations. 
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If  Mc  is  the  number  of  the  off-diagonal  terms  to  be  retained  in  a RMM  method,  and 
M x M is  the  size  of  the  full  matrix,  the  number  of  the  mathematical  operations 
needed  to  perform  a single  RMM  step  is  typically  of  the  order  of  ~ M x Mc2 
as  opposed  to  the  M3  mathematical  operations  needed  in  a single  NMM  step. 
Therefore,  the  computational  efforts  for  the  RMM  method  do  not  scale  quadratically 
with  the  increase  in  the  range  of  the  potential.  We  must  mention  that  the  RMM 
method  for  the  primary  region  dynamics,  which  we  have  used  in  this  work,  is 
independent  of  any  form  of  the  potential  in  the  primary  region  dynamics;  therefore, 
it  is  capable  of  describing  important  quantum  effects  such  as  tunneling  and  barrier 
crossing  in  the  primary  region  dynamics. 

We  have  used  a transformation  from  real  to  the  complex  time  of  the  TCFs, 
so  that  the  real  part  of  the  complex  time  always  remains  constant.  We  note  that, 
for  such  a transformation,  the  contours  of  the  complex  time  path-integral  for  the 
primary  region  dynamics,  and  that  of  the  Fourier  transform  integral  of  the  complex 
time  TCFs,  do  not  agree.  Consequently,  in  this  work,  we  would  have  needed 
information  on  many  more  points  on  the  complex  time  grid,  as  compared  to  the 
situation,  where  these  two  contours  would  overlap.  By  deforming  the  contour 
of  the  Fourier  transform  path  in  the  complex  time  plane,  we  have  derived  another 
transformation  in  which  the  contour  of  the  Fourier  transform  integral  and  the  contour 
of  the  numerical  path  integral  method  do  agree. 

We  discuss  in  what  follows  our  applications  of  our  approach,  and  suggestions 


for  further  research. 
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a)  Electronic  Absorption  Spectra  of  a Model  Two-Potential  Problem 

The  complex  time  computations  of  the  TCF  in  a single  spatial  dimension 
were  performed  by  using  the  complex  time  NMM  method.  The  results  of  this 
computation  on  one  hand  provided  a check  of  the  accuracy  of  the  numerical  method 
for  the  primary  region  dynamics,  to  be  attempted  in  the  next  application,  and  also 
a check  of  the  convergence  of  the  results  of  the  RMM  method  to  the  results  of 
the  full  NMM  method.  The  physical  and  the  computational  parameters  needed  for 
modeling  this  problem  were  the  same  as  those  of  Coalson  [Coalson,  85].  It  was 
found  that  the  real  and  the  imaginary  parts  of  the  TCF,  and  the  resulting  line  shape 
function,  computed  from  their  Fourier  transforms,  agree  very  well  with  the  results 
of  Coalson.  The  real  and  the  imaginary  parts  of  the  TCF  obtained  by  using  the 
RMM  method,  for  which  the  overall  computation  time  in  the  RMM  method  was 
reduced  by  a factor  12  as  compared  to  the  time  used  in  the  full  NMM  method, 
remained  converged  to  within  10-6  of  those  of  the  full  NMM  method.  The  only 
instructive  aspect  of  the  dynamics  that  can  be  inferred  from  these  results,  is  noted 
by  observing  a correlation  between  the  recurrences  in  the  intermediate  and  the 
long  time  behaviour  of  the  TCF,  and  the  appearance  of  the  fine  structures  in  the 
line  shape  functions.  The  structureless  smooth  line  shape  function  is  described 
primarily  by  the  short  time  behaviour  of  the  TCF,  i.e  the  Fourier  transform  integral 
for  this  type  of  line  shape  function  is  obtained  by  truncating  the  TCF  after  its 
first  decay.  The  highly  oscillating  line  shape  function  is  obtained  by  retaining  in 
the  Fourier  transform  integral  the  first  recurrence  in  the  TCF,  while  at  the  same 
time  suppressing  any  further  recurrences  that  may  be  present  in  the  longer  time 
behaviour.  The  appearance  of  higher  order  overtones  in  the  line  shape  function  is 
a confirmation  of  the  presence  of  these  further  recurrences  in  the  TCF. 
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b)  Photodissociation  of  CH3I  from  its  Ground  Vibrational  State 

Within  a TDSCF  approximation,  we  have  factored  the  TCF  for  a realistic 
polyatomic  system  (CH3I)  into  self-consistent  primary  and  secondary  region  TCFs. 
We  have  used  the  efficient  RMM  method  for  the  primary  region  dynamics,  and  have 
constructed  the  propagator  for  the  secondary  region  dynamics  analytically.  The  real 
time  behaviour  of  the  primary  and  secondary  region  amplitudes  provide  a useful 
insight  into  the  complete  dynamical  evolution  of  the  system.  At  t = 0 both  of  these 
amplitudes  start  out  as  simple  gaussian  function;  as  the  positive  real  time,  t > 0, 
increases  the  amplitudes  start  to  move  from  their  initial  positions  and  also  begin  to 
change  their  shapes.  The  dynamics  along  the  normal  mode  coordinate  X,  which 
is  dominated  by  the  dissociation  variable  R,  is  characterised  by  the  fast  oscillating 
motions  of  the  primary  region  amplitude  as  it  moves  away  from  its  initial  position. 
The  dynamics  along  the  normal  mode  coordinate  Y , which  is  dominated  by  the 
umbrella  mode  motion  variable  r,  is  characterised  by  very  slow  oscillatory  motions 
of  the  secondary  region  amplitude  as  it  also  moves  away  from  its  initial  position. 
The  contrast  between  the  fast  oscillating  motions  of  the  primary  region  amplitude 
and  the  slow  oscillating  motions  of  the  secondary  region  amplitude  are  reflected 
directly  in  the  real  time  dependence  of  the  primary  and  the  secondary  region  TCFs. 
Comparing  the  real  time  behaviour  of  the  total  TCF  with  those  of  the  primary  and 
the  secondary  region  TCFs,  we  note  that  the  real  time-dependence  in  the  total  TCF 
is  described  predominantly  by  the  real  time  behaviour  of  the  primary  region  TCF. 
The  inclusion  of  the  secondary  region  TCF  and  a time-dependent  phase  factor,  as 
described  in  chapter  4,  is  required  only  to  provide  a modulation  of  the  overall 
time-dependence  so  that  the  resultant  Fourier  transform  of  the  total  TCF  produces 
the  correct  cross-section. 
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The  theoretical  results  for  the  total  cross-section  of  the  photodissociation  of 
CH3I  from  its  ground  vibrational  state  are  compared  with  those  of  Lee  and  Heller. 
The  results  of  Lee  and  Heller  [Lee  and  Heller,  82],  which  are  obtained  by  using 
the  same  potential  energy  surfaces  as  are  used  in  our  work,  are  also  verified  by 
other  quantum  mechanical  and  semi-classical  computations  [Gray  and  Child,  84; 
Henrikson,  85;  Stodden  and  Micha,  87].  Therefore,  we  can  assume  that  for  the 
given  potential  energy  surfaces  the  results  of  Lee  and  Heller  are  reliable.  We  have 
found  a close  agreement  of  our  total  cross-section  with  those  of  Lee  and  Heller. 
Small  discrepancies  in  the  peak  position,  and  in  the  full  width  at  half  maximum, 
can  be  attributed  to  our  approximations:  expanding  the  excited  state  potential  in 
the  secondary  region  variable  Y to  retain  only  the  terms  linear  in  1 , and  using  the 
SCF  form  for  the  initial  nuclear  vibrational  state. 

Comparing  the  time-dependence  in  the  primary  and  the  secondary  region  TCFs, 
we  note  that  the  vibrations  in  the  primary  region  dynamics,  dominated  by  the 
motions  along  the  dissociation  variable  R,  are  much  faster  than  the  umbrella  mode 
vibrations  in  the  secondary  region  dynamics,  so  that  there  is  a natural  separation  in 
the  time  scales  of  these  two  vibrations  in  the  dissociating  CH3I.  It  is  perhaps  this 
separation  in  the  time  scales  of  these  two  vibrations  that  has  played  a significant  role 
in  a correct  factorization  of  the  full  dynamics  into  a primary  region  dynamics  along 
the  dissociation  variable  R,  and  a secondary  region  dynamics  along  the  umbrella 
mode  vibration  variable  r. 

Studying  the  detailed  dynamics  of  the  primary  and  the  secondary  region,  we 
note  that  at  t = 0 both  the  primary  and  the  secondary  region  amplitudes  start 
out  as  simple  Gaussian  functions,  however,  after  a sufficiently  long-time  evolution 
the  shapes  of  these  amplitudes  have  been  modulated  so  much  that  they  are  very 
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different  from  simple  Gaussian  functions.  Any  attempt,  therefore,  to  represent 
these  time  evolved  amplitudes  with  a single  gaussian  function  is  bound  to  give 
incorrect  results.  Some  of  the  approximate  quantum  mechanical  methods,  which  are 
based  on  the  the  Gaussian  Wavepacket  Dynamics  (GWD)  method  and  use  a single 
gaussian  wavepacket  for  the  dynamics,  have  also  reached  the  similar  conclusions 
[Thirumalai,  Bruskin  and  Berne,  85].  This  problem  is  generally  remedied  by  instead 
using  a collection  of  Gaussian  wavepackets,  such  that  each  of  these  wavepackets 
evolves  under  a different  set  of  equations.  However,  for  the  long  time  dynamics 
of  an  extended  molecular  system,  one  may  then  need  to  follow  the  simultaneous 
evolution  of  such  a large  number  of  coupled  Gaussian  wavepackets  that  the  whole 
computation  may  become  impractical. 

c)  Photodissociation  of  CH3I  from  Initially  Excited  Vibrational  States 

The  computations  of  the  cross-sections  for  the  photodissociation  dynamics  of 
vibrationally  excited  CH3I  were  attempted  next.  We  note  that  the  extension  of 
the  method  to  the  complex  time  evolution  of  the  primary  and  secondary  region 
amplitudes,  which  are  not  simple  Gaussian  functions,  is  straightforward.  Comparing 
the  evolution  of  the  primary  and  the  secondary  region  amplitudes,  for  the  initially 
excited  (1,0)  and  (0,1)  states  of  CH3I,  we  again  notice  a contrast  between  the 
fast  oscillations  of  the  primary  region  amplitude  and  the  slow  oscillations  of  the 
secondary  region  amplitude.  The  nodal  structure  of  the  primary  region  amplitude 
in  the  (1,0)  case  is  reflected  in  the  modulation  of  the  fast  oscillations  in  the 
corresponding  total  TCF,  whereas  the  nodal  structure  of  the  secondary  region 
amplitude  in  (0, 1)  case  causes  a modulation  only  in  the  slow  decay  of  the  total  TCF. 
The  modulation  of  the  fast  oscillations  in  the  total  TCF  of  (1,0)  case  are  reflected 
in  the  appearance  of  an  extra  peak  in  the  cross-section,  whereas  the  modulations  of 
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the  slow  decay  of  the  total  TCF  in  (0, 1)  case  appear  merely  as  a shoulder  structure 
in  the  cross-section.  This  is  a further  confirmation  of  our  earlier  stated  conjecture 
that  the  shape  and  the  position  of  the  peak  in  the  total  cross-section  is  determined 
primarily  by  the  primary  region  dynamics,  whereas  the  inclusion  of  the  secondary 
region  dynamics  causes  only  a modulation  in  the  shape  of  the  peak  so  as  to  provide 
an  overall  correct  cross-section. 

At  present,  there  are  no  fully  quantum  mechanical  computations  of  the  photodis- 
sociation dynamics  of  CH3I  from  initially  excited  vibrational  states.  The  general 
features  of  our  cross-sections  in  the  vibrationally  excited  (1,0)  and  (0,1)  cases, 
are  reported  also  by  two  other  semi-classical  computations  [Sundeberg  et  al.,  86; 
Swaminathan,  Stodden  and  Micha,  88].  However,  we  further  note  that,  while  both 
the  semi-classical  computations  have  correctly  reported  two  asymmetric  peaks  in 
the  cross-section  for  the  (1,0)  case  and  a single  asymmetric  peak  in  the  cross- 
section  for  (0, 1)  , we  could  provide  further  insight  on  the  contrasting  nature  of  the 
dynamics  along  the  dissociation  variable  R and  the  umbrella  mode  motion  variable 
r by  following  their  time-evolution. 

A detailed  study  of  the  photodissociation  dynamics  of  CH3I  from  the  other 
vibrationally  excited  initial  states  was  attempted  next.  We  fixed  the  number  of 
vibrational  quantum  in  the  secondary  region  dynamics,  and  studied  the  effect  of 
increasing  the  vibrational  excitation  energy  in  the  primary  region  dynamics.  We 
have  found  that  for  each  additional  vibrational  excitation  in  the  primary  region 
dynamics,  there  is  an  additional  peak  in  the  cross  section.  The  heights  of  these 
additional  peaks  are  not  uniform;  the  heights  of  the  peaks  in  the  high  energy 
region  of  the  cross  section,  determined  primarily  by  the  short-time  behaviour  of  the 
TCF,  are  generally  larger  than  the  heights  of  the  peaks  in  the  low  energy  region, 
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determined  mainly  by  the  intermediate  and  the  long  time  behaviour  of  the  TCF. 
This  confirms  that  the  total  cross-section  of  the  photodissociation  of  CH3I,  which 
is  a direct  process  on  a smoothly  repulsive  potential  energy  surface,  is  determined 
primarily  by  the  short  time  dynamics  in  the  Frank  Condon  region. 

As  we  put  more  and  more  vibrational  excitation  energy  in  the  primary  region 
dynamics,  we  note  a modulation  in  the  height  of  the  peak  corresponding  to  the 
secondary  region  dynamics.  This  means,  that  some  of  the  extra  vibrational  energy 
given  to  the  primary  region  dynamics  has  also  been  transferred  to  the  secondary 
region  dynamics.  A reverse  of  this  vibrational  energy  exchange,  from  the  primary 
region  dynamics  to  the  secondary  region  dynamics,  was  also  observed  by  fixing 
the  vibrational  quantum  excitations  in  the  primary  region  dynamics  and  varying  the 
excitation  energy  in  the  secondary  region  dynamics. 

The  cross-sections  for  the  photodissociation  of  CH3I  from  the  other  vibra- 
tionally  excited  initial  states,  have  not  been  verified  by  any  other  theoretical  com- 
putation or  experimental  observation.  However,  since  the  self-consistent  dynamics 
in  our  TDSCF  approximation  method  does  not  depend  upon  the  shapes  of  the 
primary  and  the  secondary  region  amplitudes,  we  expect  that  the  cross-sections 
for  these  other  vibrationally  excited  states  of  CH3I,  as  reported  in  this  work,  are 
generally  correct,  and  confirm  the  allowance  of  a self-consistent  energy  exchange 
between  the  primary  and  the  secondary  region  dynamics. 

d)  Suggestions  For  Further  Research 

Some  extensions  and  applications  of  the  TDSCF  approximation  method,  as 
developed  and  presented  in  this  work,  are  straightforward.  Simple  examples  along 
these  lines  would  be  applications  to  the  photodissociation  dynamics  of  CD3I  and 
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CF3I.  In  the  case  of  CD3I  one  can  easily  study  the  effect  of  isotopic  mass  variations 
on  the  photodissociation,  while  in  CF3I  one  has  to  incorporate  the  C-F  stretch  mode 
vibrations  besides  the  umbrella  mode  vibration  [Van  Veen  et  al.,  85].  Within  a 
TDSCF  approximation,  the  umbrella  mode  vibrations  and  the  C-F  stretch  mode 
vibrations  will  constitute  the  two-mode  dynamics  of  the  secondary  region,  whereas 
the  translational  motion  I*  along  the  dissociation  variable  will  be  treated  as  the 
primary  region  dynamics.  Since  the  experimentally  parametrized  potential  energy 
surfaces  for  both  of  these  cases  are  already  reported  in  the  literature,  the  study  of 
the  photodissociation  dynamics  of  CD3I  and  CF3I  using  the  TDSCF  approximation 
method  should  be  easy. 

Application  to  molecular  systems  in  which  there  are  more  than  one  degree  of 
freedom  in  the  primary  region  dynamics  could  also  be  contemplated.  The  most 
straightforward  example  in  this  case  would  be  the  TDSCF  treatment  of  CH3I, 
without  expanding  the  coupling  in  the  excited  state  potential  (as  was  done  in  this 
work).  The  arbitrary  general  motions  along  the  normal  mode  coordinates  X and 
Y,  would  constitute  two  degrees  of  freedom  in  the  primary  region. 

An  important  aspect  of  our  TDSCF  approach  is  the  ability  of  the  numerical 
method  to  deal  with  any  form  of  the  potential  in  the  primary  region  dynamics.  A 
possibility  would  hence  be  an  application  in  which  the  ability  of  the  method  to 
deal  with  explicit  quantum  effects,  such  as  tunneling  and  barrier  crossing  in  the 
primary  region  dynamics,  is  tested.  A recently  available  example  for  comparison, 
for  which  the  parameters  for  the  potential  energy  surfaces  and  the  results  of  a 
quantum  mechanical  close-coupling  (CC)  computation  are  listed  in  the  literature,  is 
the  photofragmentation  dynamics  of  CH3NO2  [Pemot  et  al.,  87].  In  this  test  case, 
the  dissociation  along  the  C-N  axis  in  the  excited  electronic  state  is  hindered  by  the 
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presence  of  a barrier  along  this  coordinate,  and  the  primary  region  motion  along 
the  dissociation  variable  is  coupled  to  the  harmonic  internal  motions  of  the  NO2 
fragment,  which  would  constitute  the  secondary  region  of  the  system. 

The  use  of  our  approach  for  the  long  time  evolution  of  the  amplitudes  could 
be  contemplated.  The  results  for  the  long  time  dynamics  could  be  used  to  study 
state-to-state  cross  sections  for  the  photodissociation  of  CH3I  from  its  ground  and 
excited  vibrational  states. 

The  development  of  the  formalism,  as  presented  in  this  work,  considers  dy- 
namical evolution  of  the  system  only  on  a single  potental  energy  surface,  i.e,  we 
have  not  considered  the  situation  in  which  more  than  one  potential  energy  surface 
is  responsible  for  the  dynamics  of  the  system  [Stodden  and  Micha,  87;  Coalson  88]. 
If  there  are  more  than  two  electronic  potential  energy  surfaces  in  the  formalism,  so 
that  the  evolution  of  the  amplitudes  on  any  one  of  these  potential  energy  surfaces 
is  also  coupled  to  the  other  by  the  presence  of  a non-radiative  potential  coupling 
between  the  two  surfaces,  the  frequency  operator  for  the  dynamics  of  this  coupled 
excited  surface  problem  is  written  in  a matrix  form  as 

_ [ ^11  912  1 

\Cl21  CI22 ) (8.1) 

where  flu  and  CI22  are  the  uncoupled  frequency  operators  for  the  uncoupled 

A A 1 ^ 

dynamics  on  the  two  potential  energy  surfaces,  and  ^12  = *s  the  coupling 
operator  between  the  two.  A zeroth  order  evolution  of  the  amplitudes  on  the  two 
potential  energy  surfaces,  by  ignoring  the  coupling  between  the  two,  can  still  be 
performed  through  our  TDSCF  approximation  method.  A perturbative  inclusion  of 
the  coupling  frequency  operators  f2i2  and  in  the  formalism,  causing  a mixing 


151 


of  amplitudes  during  the  dynamics,  could  also  be  attempted  in  future  work. 

Another  interesting  area  of  extension  of  the  formalism  would  be  the  inclusion 
of  anharmonic  degrees  of  freedom  in  the  secondary  region  dynamics.  To  our 
knowledge,  none  of  the  mean  field  methods  developed  so  far  try  to  include 
anharmonicity  in  the  normal  mode  vibrations  of  the  secondary  region  dynamics. 
Our  treatment  could  be  related  to  previous  work  [Vilallonga  and  Micha,  87]  on 
effects  of  anharmonicity  on  atom-polyatomic  energy  transfer.  This  would  provide 
additional  flexibility  in  the  definition  of  the  primary  and  secondary  regions. 

Lastly,  we  must  note  that  in  this  work  we  have  mainly  focussed  on  the  dynamics 
of  polyatomic  systems,  therefore,  we  have  used  the  TDSCF  approximation  for 
factoring  the  normal  mode  space  of  the  system  into  primary  and  secondary  regions. 
Since  the  general  approach,  as  presented  in  this  work,  is  not  restricted  only  to  the 
normal  mode  motions,  a possibility  would  be  to  apply  the  method  to  the  dynamics 
in  condensed  media,  where  the  dynamics  of  many  interacting  quasi  particles  could 
be  separated  into  primary  and  secondary  regions,  and  a numerical  method  could  be 
developed  for  the  exact  dynamics  in  the  primary  region. 


APPENDIX  A 

SECOND  ORDER  TRANSITION  RATES 


The  radiation  field  polarizability  tensor  TCF  for  all  second  order  processes  is 
written  as 


f*CvC  4>vc(u’0))) R • (A1) 

The  electric  field  polarizability  tensor  operator  ^(u,  0),  as  given  in  Eq.(2.15),  can 
be  expanded  as 


4>vC  (M)  = <^c+  («,  0)  + («,  0)  + *"c+  (u,  0)  + (u,  0)  (A. 2) 

where  for  example 


(u,0)  = Ev 


(+)  ^ £,(-) 


M2R  - h2u 2 


(A.3) 


the  operator  Mr  is  as  defined  in  chapter  2,  and  we  have  used  Ev  — E^  + E^~\ 
Using  Eq.(A.2)  in  Eq.(A.l)  the  polarizability  tensor  TCF,  in  the  photon  number 
representation,  is  written  as 


/&*  (“>“'»*)  = ((#C'  («'.*)*  4?  M 
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(A.4) 
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The  first  term  in  Eq.(A.4)  correspond  to  two  photon  absorption,  the  second  to  the 
fluorescent  processes,  the  third  to  elastic  and  Raman  scatterings,  and  the  fourth  is 
for  spontaneous  and  stimulated  emission  of  two  photons.  Other  cross-terms  do  not 
contribute  because  they  do  not  conserve  the  number  of  photons. 

For  example,  we  can  use  the  third  term  in  Eq.(A.4)  to  write  the  polarization 
tensor  TCF  of  the  light  source,  for  Raman  scattering  processes,  as 


Using  Eqs.(2. 18-2.20)  it  is  not  difficult  to  expand  the  TCF  in  Eq.(A.5)  as 


energy  hu £ . The  transition  rate  for  a general  Raman  scattering  process  in  the 
time  domain  description,  obtained  by  using  Eq.(A.6)  in  Eq.(2.27)  and  analytically 
evaluating  the  integrals  over  u and  u',  is  written  as 


where  Ws  is  the  statistical  distribution  of  Ns  photons  in  a mode  s = (ks,i^a)  with 


J dit+i ((x,c.(()tX,c.(0)))M  (A7) 


— oc 
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where 


(0)  ~ 


EM  + hu>r  + ie  — Hm 

k3 


-D 


C. 


(A. 8) 


is  the  molecular  polarizability  tensor  operator  for  the  target  molecule,  and  the  time 
dependence  in  this  operator  is  solely  governed  by  the  Hamiltonian  of  the  target,  i.e 


X„c.  (0  = e-*A"'X,c.  (0)  e 


+iHMt 


(A9) 


APPENDIX  B 

NUMERICAL  ANALYTIC  CONTINUATION 


The  problem  considered  here  is  the  numerical  analytic  continuation  of  a function 
f(x),  given  the  values  of  the  function  at  the  K points  x{(i  = 1,2. ..A').  One  well 
known  technique  of  numerical  analytic  continuation  by  rational  fractions  is  the  Pade 
approximant  method,  the  use  of  which,  however,  requires  the  coefficients  of  the 
Taylor  series  expansion  of  f(x).  The  Schlessinger’s  point  method  [Schlessinger, 
68],  similar  in  spirit  to  the  Pade  approximant  method,  is  instead  based  on  the 
representation  of  f(x)  by  rational  fractions,  and  uses  only  the  values  of  the  function 
at  a given  set  of  points. 

The  point  method  is  based  on  the  representation  of  the  function  at  a given  set 
of  points  by 


where 


/(*•) 


Pn  ( xj ) 
Q M (®i) 


* = I....N  + M + 1 (B.  1) 


N 

Pn(x)  = 

*=0 

M 

Qm  (ar)  = ^2qkXk  (B. 2) 

Jt=o 

Rearranging  Eq.(B.l)  one  obtains  a set  of  AT  + A/+1  linear  equations  for  jV+A/  + l 
unknown  p’s  and  q’s. 
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A more  efficient  way  to  effect  the  point  type  of  representation  is  obtained  by 
using  continued  fractions.  To  do  this,  consider  the  continued  fraction 


Cn  (x)  = 


/ 00 


1 + ai  (x  - x i ) 


1 + a2  (x  — X2 ) 


1 + 


(5.3) 


QN  (x  ~ Xtf) 


It  is  easy  to  see  that  Cn(x\)  = /(x i)  and  that  the  coefficients  a\...as  can  be 
chosen  so  that  C/v(x,)  = /(x,).  The  coefficients  aj  in  Eq.(B.3)  are  determined 
recursively  from  the  formula 


a 1 { 1 + aj—i(xj+ 1 xj—i) 

( Xj  xj+ 1)  1 + dj  — 2 (xj^-l  Xj— 2) 


1+  ... 


(5.4) 


, aj(xj+ 1 xj ) -i 

•"+  /(x.)  * 

1 TfcTTT 


with 


{[f  M / f M]  - 1} 

(x2  - Xi) 


a 1 = 


(5.5) 


APPENDIX  C 

SCF  EIGENVALUES  and  EIGENFUNCTIONS 


SCr  coefficients 
CP 

0.99672584059640 
8 . 0792096484  369D-02 
3 . 1993488482456D-03 
EP 

5 . 267 51 092854 49D-03 
E(0,0) 

SCr  coefficients 
CP 

8 . 20391 0854 3589D-02 
-0.97791519095421 
-0.19222763581328 

EP 

1.0195151579200D-02 
E( 1 , 0 ) 

SCF  coefficients 
CP 

0.99294262572004 
0.11822018162167 
9 .4303067565885D-03 

EP 

1 . 08602455237B5D— 02 
E(0,1) 

SCF  Coefficients 
CP 

0.11876956673886  -2, 

-0.96472325967676  -0 

-0.23495280857082  4 

EP 

1.57  46 5764886 BID-02  2 
E(l,l)  - 2 

SCF  Coefficients 
CP 

1 . 579B871824954D-02 
-2 . 4974466984030D-02 
0.32261114169622 
-0.84782678341756 
-0.41980765622034 
EP 

1 .99984 31 071516D-02 
E( 3 , 0 ) 

SCF  Coefficients 
CP 

-1.1413194075059D-02 
0.25203921154663 
-0.88546300889190 
-0.38277108867786 
-7 . 6075806932522D-02 
EP 

2 . 059 50412866 36D-02 
E(2,l) 

SCF  Coefficients 
CP 

-5 . 1220220292604D— 03 
0.20241294293276 
-0.91948962493189 
-0.33269077247215 
-5.34644272215150-02  - 

EP 

1.50731147119770-02 

1(2.0)  - 


for  (0,0) 

CS 

0.99992129098795 
-1.23386486956770-02 
-2 . 2736704580636D-03 
ES 

6 . 7930625343156D-03 
8.0436394067040D-03 
for  (1,0) 

CS 

0.99975548326356 
■2. 20851859133370-02 
-1.1037425480124D-03 
ES 

9.0966254686394D-03 
1. 296 396 139886 3D-02 
for  (0,1) 

CS 

-1.5912198657  396D— 02 
-0.99944224407830 
2.93598822479530-02 
ES 

1.7891976191543D-02 
1.91397468013950-02 
for  (1,1) 

CS 

. 3940237480153D-02 
.99888792882731 
. 0617  36909B342D-02 
ES 

. 00954 2175791 0D-02 
.4032598177014D-02 
for  (3,0) 

CS 

-0.99967786508303 

2.53341804931680-02 

9.2643352437270D-04 

1.2176546302277D-03 

-6.6331753138563O-05 

ES 

1.2851509975508D-02 
2 . 2766312540772D-02 
for  (2,1) 

CS 

-3 . 5129412140848D-02 
-0.99779421024850 
5.6238813365018D-02 
-1.7967262680218D-03 
2.5702224089621D-03 
ES 

2.2432566310398D-02 
2.8B66752078828D-02 
for  (2,0) 

CS 

0.99946258622104 
3.2753228413340D-02 
3 . 99482 51 3300 55D-04 
1.2641843950966D-03 
8.3811681449883O-05 
ES 

1.14475261609390-02 

1.7832885380423O-02 


EPS 

-4.0169340561566D-03 


EPS 

—6 . 3278 156 4 89 76 2D— 03 


EPS 

-9 . 6124749139322D-03 


EPS 

1.1809400069577D-02 


EPS 

-1.00836285062S1D-02 


EPS 

-1 .41608555182070-02 


EPS 

-B.6877554924933O-03 
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SCF  Coefficients 
CP 

1 . 639066 198 3180D— 02 
-3.93997969963300-02 

0.37176544438281 

-0.81036971037964 

-0.45083868294623 

EP 

2 . 5 S3 08 02 2191 88D— 02 
E ( 3 , 1 ) 

scr  coefficients 
CP 

-0.98732804812872 
-0.15586673870250 
-2. 19481384417870-02 
-1.9366257119792D-02 
-5.6667846867243D-03 
EP 

1 . 6366745947212D— 02 
E ( 0 , 2 ) 

BCF  Coefficients 
CP 

-0.98043470108244 

-0.19275011449047 

-3.2853518085547D-02 

-2.1557264133552O-02 

-7.1498943887409D-03 

EP 

2 . 1922 4993 52493D-02 
E( 0 , 3 ) 

SCr  coefficients 
CP 

0.15908194189023 
-0.93866044635406 
-0.29658955223327 
-6 . 3944553417524D-02 
-3.9436452795984O-02 
EP 

2.1239326215826D-02 
E(l, 2)  - 

scr  coefficients 
CP 

0.19493661408523 
-0.91499001974979 
-0.34015682865364 
-8 . 4662836984030D-02 
-4. 3799979756808O-02 
EP 

2 .66 4280609753 3D-02 
E(l, 3)  - 

SCr  Coefficients 
CP 

6 .4410500860338 D-03 
-3 . 6266  39209693BD-02 

8.11509506542620-02 

-0.39563982745886 

0.91407161197266 

EP 

2.6825547377424O-02 

8(6,0) 


for  (3,1) 

cs 

-2.1545029807526D-02 
-0.99906974566358 
3.7208893384953D-02 
2 . 290589981 356 4D-03 
2 . 368815795899 5D-03 
ES 

2.3743167197191D-02 
3 . 384 39990794 46D-02 
for  (0,2) 

CS 

-1.80586610518990-03 
-3 . 4519753651433D-02 
-0.99811648229037 
5.0628486759489D-02 
2 . 3172320095829D-03 
ES 

2.8999861B18118D-02 
3 . 021 55711S9619D-02 
for  (0,3) 

CS 

-1.336  392120206SD-03 
-2.8461322257134D-03 
-5.66857819881730-02 
-0.99552466063534 
7.5397423443674O-02 
ES 

4. 00749024042770-02 
4 .12814564460890-02 
for  (1,2) 

CS 

■1 .14  385267092730-03 
■5. 06037701021810-02 
■0.99625581324658 
7 . 0067 4090 31740D-02 
•1.6919858235142  D-03 
ES 

142298 08674 15D-02 
3.5046603253163D-02 
for  (1,3) 

CS 

1.2880436744038D-03 
•1.9539643892144 D-03 
•7. 5989741313275O-02 
■0.99243908846012 
9.6357344703184D-02 
ES 

4 .246504 962427 3D-02 
4.60775225110440-02 
for  (4,0) 

CS 

-0.99790811482856 
-6.3779166082072D-02 
1 . 0233873792702D-02 
2.6226699935523O-03 
4 .19818671029660-05 
ES 

1.9641S98824394D— 02 
2.96886526673390-02 


EPS 

-1 . 5429970336933D-02 


EPS 

-1.5151036605712D-02 


EPS 

-2.0715945310681D-02 


EPS 

1.7615703830078D-02 


EPS 

2 . 3030333210761O-02 


EPS 

-1.6778493534479O-02 
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